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Abstract: Dark matter in the sub-GeV mass range is a theoretically motivated but largely unex¬ 
plored paradigm. Such light masses are out of reach for conventional nuclear recoil direct detection 
experiments, but may be detected through the small ionization signals caused by dark matter-electron 
scattering. Semiconductors are well-studied and are particularly promising target materials because 
their 0(1 eV) band gaps allow for ionization signals from dark matter particles as light as a few hun¬ 
dred keV. Current direct detection technologies are being adapted for dark matter-electron scattering. 
In this paper, we provide the theoretical calculations for dark matter-electron scattering rate in semi¬ 
conductors, overcoming several complications that stem from the many-body nature of the problem. 
We use density functional theory to numerically calculate the rates for dark matter-electron scattering 
in silicon and germanium, and estimate the sensitivity for upcoming experiments such as DAMIC and 
SuperCDMS. We find that the reach for these upcoming experiments has the potential to be orders 
of magnitude beyond current direct detection constraints and that sub-GeV dark matter has a sizable 
modulation signal. We also give the first direct detection limits on sub-GeV dark matter from its scat¬ 
tering off electrons in a semiconductor target (silicon) based on published results from DAMIC. We 
make available publicly our code, QEdark, with which we calculate our results. Our results can be 
used by experimental collaborations to calculate their own sensitivities based on their specific sefup. 
The searches we propose will probe vast new regions of unexplored dark matter model and parameter 
space. 
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1 Introduction 

1.1 The search for dark matter 

There has been tremendous progress in the last three decades in the direct detection search for weak- 
scale dark matter (DM) using underground detectors. The original aim was to probe the scattering 
through Z-exchange of DM candidates with roughly weak-scale mass against nuclei [1]. Now, exper¬ 
iments searching for these DM-induced nuclear recoils [2-4] are sensitive to scattering cross sections 
many orders of magnitude below the Z-exchange cross section, for candidates in the O(10GeV- 
lOTeV) mass range. The motivation behind this incredible experimental achievement has been the 
theoretically appealing, and dominant. Weakly Interacting Massive Particle (WIMP) paradigm: DM 
as a weak-scale thermal relic associated with new physics that solves the hierarchy problem. However, 
the era of this paradigm’s preeminence appears to be ending due to both the lack of a DM discovery, 
which excludes significant regions of WIMP parameter space, and the absence of non-Standard Model 
(SM) physics at colliders, which has undermined the theoretical motivation behind it. More impor¬ 
tantly, several other theoretically motivated candidates exist for resolving this great mystery of particle 
physics. 

Motivated particle-DM candidates have been proposed over a vast range of masses, from ultra¬ 
light bosonic fields such as a QCD axion [5-7], to non-thermal GUT-scale relics [8]. While these have 
inspired a diverse array of experimental searches, techniques for probing them are far less developed 
than the WIMP search program. One well-motivated candidate that has received increased attention 
recently and is the focus of this paper is light dark matter (EDM), with DM masses in the MeV to GeV 
range. EDM is often motivated by production mechanisms that go beyond the standard freeze-out and 
may be found in several frameworks in which the sub-GeV mass scale arises naturally. In addition, 
the origin for the DM relic density can be naturally addressed by several mechanisms that suggest that 
EDM interacts with SM particles via, for example, an exchange of a light “dark photon”, an axion, or 
through an electromagnetic dipole moment. There is a large range of parameter space of such models 
that evades both laboratory and astrophysical bounds [9-30]. 
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Investigating LDM is an important and natural direction to pursue in the DM search effort. An 
essential part of this pursuit is extending direct detection searches to this low mass range. Several 
possible ways to do this were described in [9]. Fortunately, much of the impressive technology being 
developed for the Weak-scale direct detection program can be readily adapted to search for LDM. An 
example of this was described in [31], obtaining the first direct detection limits on DM with masses 
as low as a few MeV using published XENONIO data. In this work, we study in detail the even more 
promising possibility of semiconductor-based LDM searches, significantly expanding the preliminary 
work done in [9]. Other, complementary techniques to search for LDM have been discussed in [32- 
61]. 

1.2 Direct detection of sub-GeV dark matter 

Current direct detection experiments are limited to probing DM masses above a few GeV due to the 
high energy thresholds required for detecting nuclear recoils. The challenge in probing lower DM 
masses is twofold: for lower masses, not only is the total kinetic energy of the DM particle decreased, 
but so is the fraction of energy that is transferred to the nucleus. As a result, the energy of the nuclear 
recoils is much lower and one must drastically reduce the threshold energies to detect it. This is an 
experimentally challenging task, although it may be possible to probe masses down to a few hundred 
MeV, see [62, 63]. Instead, as discussed in [9], scattering channels other than elastic nuclear recoil 
are likely to be far more fruitful. 

A very promising avenue is to search for the small ionization signals caused directly by DM- 
electron scattering. The lightness of the electron and the inelastic nature of the DM-electron scattering 
process allow DM particles to transfer a large fraction of their kinetic energy to the electron when 
they scatter, enabling DM as light as ~1 MeV to cause an ionization signal. Furthermore, detecting 
small ionization signals is already a well-developed part of direct detection technology. In fact, the 
XENONIO experiment was already sensitive to the ionization of a single electron [64], and results 
of a short single-electron-sensitive run [65] were used in [31] to place direct detection bounds on 
DM with masses as low as a few MeV. This serves as a proof-of-principle, motivating dedicated 
EDM searches in other dual-phase noble liquid experiments such as XENONIOO and EUX. However, 
semiconductor targets have the potential to probe even smaller cross sections. In semiconductors such 
as silicon or germanium, the band gap (the threshold to “ionize” an electron by exciting it from a 
valence band to a conduction band) is ~1 eV - a factor of 10 to 20 times lower than the ionization 
threshold in liquid xenon. The consequences of this lower energy threshold are significant. Not 
only could this allow sensitivity to DM down to masses below an MeV, but it would also mean a 
substantial increase in event rate for all DM masses [9, 25]. The reason for this is that, given the 
characteristic velocities of DM particles and electrons, ~1 eV recoil energies are typical, while recoil 
energies of ~10eV require velocities that are only found on the tails of the DM and electron velocity 
distributions. Moreover, although the background that causes an ionization signal at such low energies 
is still poorly understood, it is reasonable to expect that background event rates in semiconductors 
may be significantly lower than in xenon-based detectors [66] (especially since they may be operated 
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cryogenically)*. There is currently an active program in, for example, both the SuperCDMS and 
DAMIC collaborations to develop germanium- and silicon-based detectors that are sensitive to single 
electron-hole pairs [66, 67], enabling a leap forward in LDM detection. 

This developing experimental program presents a new theoretical challenge: the calculation of 
the expected signal rate. Unlike for elastic nuclear recoils, this calculation is highly non-trivial. In this 
paper, we tackle this calculation head on and present detailed new results for germanium and silicon 
targets. 

1.3 The challenge of calculating event rates 

Several factors complicate the calculation of DM-electron scattering rates. Bound electrons in dense 
media have a) typical speeds of order a k, 1/137 or greater, much faster than DM particles (with 
V ~ 10“^), b) indefinite momentum, with even very large momenta having non-zero probability, and 
c) a complicated structure of energy-levels. This greatly modifies the scattering kinematics and breaks 
the simple link between momentum transfer and energy deposition. As we discuss in more detail 
below, event rates can be highly sensitive to the energy-level structure and the tails of the electrons’ 
momentum distributions. In addition, the quantum nature of both the initial and final elecfron stales 
is important, and they cannot be correctly treated classically. As a result, approximate calculations 
which do not fully account for these details may not give accurate results. This becomes even more 
important for the large energy depositions, well above OfeV), since these rely on the tails of the 
electron’s momentum distribution. Once correctly calculated, the effect of all these complications 
can be completely encoded in an atomic form factor [9]. This function is different for each specific 
target material, but is independent of the DM model. Once it is known, event rates can be calculated 
relatively simply. 

When the target is an isolated noble gas atom, the combination of spherical symmetry and 
previously-compiled bound-state electron wavefunctions makes calculation of the ionization form 
factor relatively straightforward. Refs [9, 31] used this as an approximation for the form factor of a 
liquid xenon target. However, calculating the form factor for a crystal target (such as a semiconductor) 
is far more challenging. A periodic crystal lattice is a complex multi-body system, with outer-shell 
(valence) electrons delocalized and occupying a complicated energy band-structure. Accurate wave- 
functions of the valence electrons cannot be found analytically, but must be computed numerically 
with an expansion in a discrete set of plane waves^. Taking this approach, a first calculation was done 
in [9], assuming a single-electron threshold in a germanium target. A second approach was taken 
in [25], which succeeded in simplifying the calculation until it was analytically tractable. However, 
the approximations required for this were so extensive that the result might be considered only as an 
order-of-magnitude estimate. A third, semi-analytic approach was taken in [68] (see “Note Added”), 
where numerical bound-state wavefunctions for free germanium and silicon atoms are used and the 

'Unlike for traditional WIMP searches, nuclear recoils are not an important background for our electron recoil signal as 
their rates are expected to be much lower than background-induced electron recoils. 

^Note that inner-shell (core) electrons, which are important in some cases, are more localized so that their wavefunctions 
are closer to those computed assuming an isolated atom. 
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Figure 1. Prospects for Upcoming DM-Electron Scat¬ 
tering Searches: Selected near-term projections of DM- 
electron scattering cross-section as a function of DM 
mass my^ for the DAMIC (green curves) and SuperCDMS- 
silicon (dark red curves) experiments, for different ioniza¬ 
tion thresholds and (hackground-free) exposures, as indi¬ 
cated. Solid curves show the 95% C.L. exclusion reach 
from simple counting searches, while dashed curves show 
the Scr-discovery reach from annual modulation searches. 
The gray shaded region shows the current XENONIO 
bound [31], while the shaded green region shows the esti¬ 
mated bound from 2012 DAMIC data with a ~11-electron- 
hole pair threshold. The projections for SuperCDMS- 
germanium (not shown) are comparable to silicon. See §6.5 
for more details. The three plots show results for the differ¬ 
ent indicated DM form factors, corresponding to different 
DM models. 


outgoing electrons are described by plane waves. The latter approach gives answers much closer to 
our full numerical calculation, but important differences remain. 

1.4 Overview of the paper 

In this paper, we present the results of a detailed numerical computation of DM-electron scattering 
rates in germanium and silicon targets as a function of the electron recoil energy. This significantly 
expands on the previous calculation in [9]. Higher recoil energies for the scattered electron allow 
a larger number of additional electron-hole pairs to be promoted via secondary scattering. Using 
a semi-empirical understanding of these secondary scattering processes, we convert our calculated 
differential event rate to an estimated event rate as a function of the number of observed electron-hole 
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pairs. These results will allow several experimental collaborations, such as DAMIC and SuperCDMS, 
to calculate their projected sensitivity to the DM-electron scattering cross-section, given their specific 
experimental setups and thresholds. It will also allow them to derive limits on this cross section in the 
absence of a signal, or the preferred cross section value should there be a signal, in forthcoming data. 
Achieving low ionization thresholds could allow these experiments to probe large regions of LDM 
parameter space in the near future, as illustrated in Fig. 1. 

In §2, we briefly discuss the direct detection prospects for a few popular LDM models. We will 
see that the upcoming generation of experiments with semiconductor targets play an essential role in 
testing these models. In §3, we outline how to calculate the rate for DM to scatter off bound electrons. 
We provide an intuitive understanding of the scattering kinematics. Our discussion is general and 
applicable to both electrons bound to (free) atoms as well as electrons in semiconductor targets. The 
details of this calculation as well as comprehensive formulas are contained in Appendix A, signifi¬ 
cantly expanding on the information contained in [9, 31]. We then focus on semiconductor targets, 
and describe the numerical computation of the scattering rates in §4. We describe our code QEdark, 
which is an additional module to the publicly available code Quantum ESPRESSO [69]. The latter 
calculates the band structure and all electron wavefunctions using density functional theory (DFT) 
and pseudopotentials, two established condensed matter computational tools, to calculate the Bloch 
wavefunction coefficients for the initial and final state electrons. In QEdark, we use this information 
to calculate the crystal form factor for DM-electron scattering as well as the scattering rates. QEdark 
and the crystal form factors will be publicly available at this link. In §5, we discuss the conversion 
from the energy of the primary scattered electron to the size of the final ionization signal. We present 
a conversion formula and discuss the uncertainty associated with it. In §6, we present the results of 
our computation, showing the cross-section sensitivity as a function of detector threshold, as well as 
the potential discovery reach using annual modulation. We also provide detailed sensitivity estimates 
for two representative, near-term experiments that may soon reach the required sensitivity to detect 
LDM, namely DAMIC and SuperCDMS. We conclude in §7. The appendices contain additional tech¬ 
nical details: Appendix A provides a detailed derivation of the formulae for the scattering rate and 
crystal form-factor. Appendix B describes our choice of local DM velocity distribution. Appendix C 
discusses the convergence of our numerical results. Appendix D studies the effects of inner-shell elec¬ 
trons on the overall scattering rate. Appendix E presents details of the systematic study of secondary 
interactions, and Appendix F gives a brief review of DFT and pseudopotentials. 

We note that our main results are contained in Figs. 1, 2, 6, and 9 and described in §2 and §6. 

2 Models of Light Dark Matter 

Theories of LDM have been receiving increased attention in recent years. Here we illustrate with 
just a few benchmark LDM models how the upcoming generation of experiments with semiconductor 
targets, including SuperCDMS and DAMIC, play an essential role in the search for LDM. Classes 
of models that are probed by LDM direct detection include DM that scatters through a dark-photon 
mediator or through a dipole moment interaction. We focus on DM coupled to a dark photon, leaving 
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a discussion of dipole moment interactions [25, 70], the SIMP [29, 30], and other models that can be 
constrained by electron recoils to an upcoming publication [71]. 

For illustration, we consider models of LDM based on the vector-portal, in which the dark sector 
(and the DM particle, x) communicates with the SM through a ?7(1 )d gauge boson A'. The A' is 
kinetically mixed with the SM hypercharge U{1)y via the interaction 


£ D 


_f_ 7?/^^ /?' 

2COS0W ^ 


(2.1) 


causing it to couple dominantly to electrically charged particles at low energies. Here e is the kinetic 
mixing parameter, dw is the Weinberg mixing angle, and i^ ih® field 

strength. 

DM particles can scatter off electrons in direct-detection experiments through A' exchange. In 
the notation of §3.2 below, the DM-electron reference cross section is given by 
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where is the DM-electron reduced mass and an = (with go the U{1)d gauge coupling). 

We note that this expression is the same for DM that is a complex scalar or a fermion. The corre¬ 
sponding DM form factor is 


Fouig) 


m\, -I- a^ml f 1, niA' > amg 
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(2.3) 


where q is the momentum transfer between the DM and electron. 

In Fig. 2, we illustrate the parameter spaces of both the uia' S> arUg and niA' arUg regimes, 
taking the fermionic and complex-scalar cases separately for the former. We study three cases, which 
highlight different possible production mechanisms, and show the interplay between different experi¬ 
mental probes. 


(i) Freeze-out via the vector portal: complex scalar LDM 

We consider the phenomenologically interesting and predictive region > 2m^, corre¬ 
sponding to Fouiq) = 1- Annihilation to SM particles occurs via an off-shell A' {xx* 

A'* -A SM). This process is p-wave suppressed, allowing the DM abundance to be set by 
thermal freeze-out while evading constraints from the cosmic microwave background (CMB), 
e.g. [72, 73], and from gamma-rays in the Milky-Way halo [74]. We show the parameter space 
for this scenario in Fig. 2 {top left), taking for concreteness. The thick blue 

curve shows the cross-section for which the correct relic abundance is obtained from freeze- 
out [73] (this is largely insensitive to the specific choice of niA')- Above this line, an asym¬ 
metric DM component may complete the DM abundance. Below it, the abundance is naively 
too large, but this region may be viable with alternate hidden-sector freeze-out channels. We 
also show various constraints on this model. The black curve labelled “XENONIO” shows the 
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Figure 2. Prospects for Benchmark Models: Selected 
95% C.L. exclusion reach for the DAMIC (green curves) 
and SuperCDMS-silicon (dark red curves) experiments, 
compared with other constraints for the benchmark mod¬ 
els discussed in §2. White regions are unconstrained, while 
thick blue curves illustrate possible predictive mechanisms 
for generating the DM abundance. Top: DM interacting 
via a massive dark photon {Fumiq) ~ 1), for complex- 
scalar DM with freeze-out abundance (left), and Dirac- 
fermion DM with asymmetric abundance (right). Bottom: 
DM interacting via an ultralight dark photon (fbM(5) = 
(ame/g)^), with an abundance generated by freeze-in. The 
DAMIC and SuperCDMS projections assume 100 g-year 
and 10 kg-years background-free exposures, with 2- and 1- 
electron thresholds, respectively, in a silicon target. See text 
for details. 


electron-recoil DM constraint set with XENON 10 data [31]. The black curve labelled “Current 
NR Constraints” shows constraints from conventional nuclear-recoil searches from [3, 75, 76]. 
Some measurements only constrain e as a function of uia'- Among these, we only show the 
strongest constraints, which are a BaBar search for e+e“ —)• 7 -F invisible [49, 51, 52] as well 
as electroweak precision tests (EWPT) [77, 78]; however, to guide the eye, we also show the 
“favored” 2 (T-region for which the A' can explain the discrepancy between the measurement 
and SM prediction for the muon anomalous magnetic moment, [79]. We translate these 
into the ag versus plane by using the constraint on a d from either perturbativity [80] or x 
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self-interactions [22]. For these we require that is less than 1.0 and small enough so that 
cJseif-int/^^x ^ 1 cm^/g for clusters [81]. A second set of constraints bound some combina¬ 
tion of e, od, and rrijii: the electron beam-dump E137 [57, 82] and the proton beam-dump 
LSND [60, 83, 84]. We again use the constraint on from self-interactions and perturba- 
tivity to translate these into the ae versus plane. We also show a rough bound on N^s, 
see [53, 54, 73]; the presence of additional relativistic degrees of freedom could allow this 
bound to be evaded. For a complementary representation of this parameter space see [59]. 

(ii) Freeze-out via the vector portal: Dirac fermion LDM 

In Fig. 2 {top right), we consider the same scenario as in (i) but take x to be a Dirac fermion. 
This also corresponds to FbM(7) = 1- The main difference between this scenario and (i) is that 
the annihilation cross section is now s-wave, so that constraints from the CMB preclude the 
abundance being set by freeze out. Instead, we assume the abundance to be asymmetric [85- 
87], and require the symmetric component to be small enough after freeze-out to avoid the 
CMB bounds [20]. This provides a lower bound on the annihilation cross-section and thus on 
(7e, shown with a black solid line. As before, this lower bound is model-dependent and can be 
evaded with additional annihilation channels. The other constraints are similar. 


(iii) Freeze-in via the vector portal 

In Fig. 2 {bottom), we consider an ultra-light A' mediator {mA’ ^ came), corresponding to 
Fdm{<i) = {aruelq)'^- Here the couplings are so small that the DM would never have thermal- 
ized with the SM sector. The x abundance can receive an irreducible “freeze-in” [88] contribu¬ 
tion from 2 —)• 2 annihilation of SM particles to xx well as Z-boson decays to XX’ computed 
in [9] (see also [24]). The parameters required for the abundance again uniquely constrain (ff, 
versus m^, as shown by the thick blue curve. In addition to the XENONIO electron-recoil con¬ 
straint [31], we also show the bounds from conventional nuclear-recoil searches. The nuclear 
recoil cross-section, (Jnr, can be related to the electron recoil cross-section by 


d{(7NRv) 

di^NR 


Z‘^ae{ame)‘^ 


7(^min,NR) ; 


(2.4) 


where the target nucleus has mass m^v and atomic number Z, T^nr is the nuclear recoil energy, 
and Umin,NR = (2/x^Ar), V is relative velocity of the DM, and rj is the inverse mean 

speed defined in Appendix B. Since this recoil spectrum is peaked towards low energies more 
than for a contact interaction, determining accurate DM constraints requires a careful analysis 
of the experimental data. We place approximate bounds from “CDMSEite” [89] and EUX [3] 
results, taking the former to have 6.2kg-days germanium exposure, a 0.84keV threshold, 100% 
signal efficiency and 10 observed events, and taking the latter to have 10 tonne-days xenon 
exposure, a 5 keV threshold, 50% signal efficiency and 0 observed events. Due to the smallness 
of the couplings, the other constraints seen in the previous scenarios are absent in this one. 
Instead, we include various astrophysical constraints on millicharged particles, which are also 
applicable for DM coupled to an ultralight A' [90]. 





In each of the figures in Fig. 2 we show the prospects for DAMIC (100 g-years, silicon target, 
2-electron threshold) and SuperCDMS (10 kg-years, silicon, 1-electron threshold), discussed in §6.5. 
We note that a magnetic-dipole-moment interaction would also give Fdm{q) = while an electric- 
dipole-moment interaction would give Fdm(<?) = amejq- We see that these models above all have 
concrete predictions that the upcoming generation of direct detection experiments can test. 


3 Direct detection of dark matter by electron scattering in semiconductors 

In this section, we review the theory of DM scattering with bound electrons. We begin in §3.1 by 
considering the simple kinematics of LDM scattering with both nucleons and electrons. This makes 
clear the motivation for using electron recoils to probe LDM. The discussion also shows that the 
DM-electron scattering rate is expected to be sensitive to the details of electron binding in the target, 
especially for higher energy/ionization thresholds. A consequence of this is that to calculate accu¬ 
rate scattering rates, detailed modeling of the electronic structure of the target material is required, 
involving knowledge of the wavefunctions of all accessible occupied and unoccupied electron levels. 

In §3.2, we summarize how this scattering-rate calculation is formulated, with a focus on the case 
of semiconductor targets. The key results are Eqs. (3.13) and (3.17). The former gives the differential 
scattering rate in terms of the DM model, the DM velocity profile, and crystal form factor. The 
latter gives the crystal form-factor, which encodes all the relevant electron binding effects for a given 
target material. This reviews and extends the discussion from Ref. [9]. In Appendix A, we provide a 
full derivation of all the results given here. For the interested reader, in Appendix A.3 we present a 
derivation of the ionization rate in free atomic targets, as is relevant for xenon targets and which was 
used in Refs. [9, 31]. 


3.1 Kinematics of dark matter-electron scattering 

Conventional DM direct detection experiments assume that the DM particle scatters elastically off a 
target nucleus. This recoiling nucleus then collides with the surrounding matter within the detector, 
giving off energy in the form of heat, phonons, ionized electrons, scintillation photons, etc, depending 
on the detector material. However, if the DM particle is light, the momentum transfer, q, between the 
DM and the target nucleus is small and may not provide enough energy for the recoil of the nucleus 
to be detected. We can see this through the following argument. The energy of the recoiling nucleus 
in nuclear scattering is 


.E'nr 


g' < ~ 1 eV X ( f . 

2mAr rriN VlOOMeV/ \ m^v / 


(3.1) 


For the scaling in the last step of this equation, we have taken the typical DM speed to be 300 km/s 
« 10“^c, and assumed <C rriN- For = 30 GeV, we find F^nr ~ 2 keV. However, if we 
consider lighter DM masses, such as = 100 MeV, the recoil energy drops to Snr ~ eV, which 
is well below the detection thresholds of current direct detection experiments (e.g. ~ 840 cVnr for 
CDMSlite [89] and ~ 4 keVNR for LUX [3]). Note that the energy of the recoiling nucleus is also not 
efficiently transferred to electrons, and so is not nearly large enough to ionize or excite even a single 
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electron; it is also well below current phonon detection thresholds. As a result, DM masses below a 
few hundred MeV escape detection no matter how large their cross section. 



Figure 3. The scattering of a DM particle with a bound electron. The DM transfers momentum q to the target, exciting it 
from the ground state X to an excited state X*, which can be either a higher-energy bound state or an ionized state. 


Now consider a DM particle colliding directly with a bound electron, exciting it to a higher 
energy level or an unbound state, as illustrated in Fig. 3. The kinematics are very different from those 
of a nuclear recoil. Firstly, being in a bound state, the electron does not have definite momentum - 
in fact it may have arbitrarily high momentum (albeit with low probability). This breaks the direct 
relation between recoil energy and momentum transfer given in Eq. (3.1). The energy transferred to 
the electron, AE^, can still be related to the momentum lost by the DM, q, via energy conservation: 


AEe = -AE^ - AEm = ^ - 

^ 2mv 2 ^ 


2m M 


= q ■ V — 




(3.2) 


Here the AE^ term accounts for the fact that the whole atom also recoils. In practice this term is 
small, which also allows us to replace with m^. We thus define 


Ee = AEe = -AE^ 


(3.3) 


as fhe energy fransferred fo fhe elecfron.^ Since an arbifrary-size momenfum fransfer is now possible, 
fhe largesf allowed energy fransfer is found by maximizing AEg wifh respecf fo q, giving 


1 n 1 

AEe < 2FxJvT ~ - eV X 


/ \ 

\MeVJ ■ 


(3.4) 


This shows fhaf all fhe kinetic energy in fhe DM-afom collision is (in principle) available fo excife fhe 
election. For a semiconducfor wifh an 0(eV) bandgap, ionizafion can be caused by DM as lighf as 
C>(MeV). 

Whaf is fhe likelihood of acfually obfaining a large enough q fo excife fhe election? This brings 
us fo fhe second major difference compared fo DM-nuclear scattering: fhe election is bofh fhe lighfesf 
and faslesl parficle in fhe problem. The fypical velocity of a bound election is Ve ~ Z^ga, where 
Zeff is 1 for oufer shell elections and larger for inner shells. This is much greafer fhan fhe fypical DM 

^We emphasize that Ee is the energy transferred to the electron, not its kinetic energy. Some of this energy goes 
to overcoming the binding energy. As we will discuss further in §5, in semiconductors the remaining energy is rapidly 
redistributed by secondary scattering processes, which can produce further electron-hole pairs. 
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velocity of ~ 10 The typical size of the momentum transfer is therefore set by the electron’s 
momentum, 

qtyp f^xe^rel - meVe ~ ^ ^eff X 4 keV (3.5) 


where Urei is the relative velocity between the DM and electron. 

Returning to Eq. (3.2), the first term on the right dominates as long as is well above the bound 
in Eq. (3.4). This gives a formula for the minimum momentum transfer required to obtain an energy 
AEe. 


q ^ 




AE, 


qtjp 


(3.6) 


4 Zes eV 

This scaling suggests that the typical available momentum is enough to cause a transition of just a 
few eV, such as for an electron being excited just across the germanium or silicon bandgap. Exciting 
a more energetic transition will require a momentum out on the tail of the electron’s momentum- 
space wavefunction (or probing the tail of the DM velocity distribution), and its probability will be 
correspondingly suppressed (as can be seen clearly in Eig. 5 below, which we will discuss in §6.1). 
Ionization of a xenon atom, requiring ~10eV energy, falls into the second category, as do most 
possible transitions to the conduction band in germanium or silicon. 

Erom this argument we expect the rate of DM-electron scattering to be sensitive to the precise 
forms of the electron energy levels and wavefunctions in the target. The computation we present below 
is designed to address this sensitivity by modeling in detail the electronic structure in germanium and 
silicon crystals. A corollary of this argument is that, given the u-dependence in Eq. (3.6), the rate 
should also be sensitive to the DM velocity profile. As this varies over the year, we expect a significant 
annual modulation in the signal size, a potentially crucial test of the DM origin of a signal. We discuss 
the expected annual modulation in §6.4. 


3.2 Calculating excitation rates 

3.2.1 General formulation for dark matter-induced electron transitions 

If a DM particle scatters with an electron in a stationary bound state such as in an atom, it can excite 
the electron from an initial energy level 1 to an excited energy level 2 by transferring energy AEi ^2 
and momentum q. The cross section for this process takes quite a different form to the free elastic 
scattering cross section. 

If Alfree(9) is the matrix element for free elastic scattering of a DM particle and an electron, then 
we parametrize the underlying DM-electron coupling using the following definitions [9] : 

\Miree{qW = |Affree(a^e)^ X |TbM(9)l^ 

_ _ /i2g|Affree(ame)P 
167rm|m^ ’ 

where | Af P is the absolute square of Af, averaged over initial and summed over final particle spins. 
The DM form factor, AdmI?). gives the momentum-transfer dependence of the interaction - for ex¬ 
ample, EDyi{q) = 1 results from a point-like interaction induced by the exchange of a heavy vector 


(3.7) 

(3.8) 
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mediator or magnetic dipole moment coupling, FomC?) = {otme/q) for an electric dipole moment 
coupling, and Fdm(9) = {ocme/qY for exchange of a massless or ultra-light vector mediator (see 
§2). Ue parameterizes the strength of the interaction, and in the case of FomI?) = 1 is equal to the 
cross section for free elastic scattering. All sensitivity estimates or constraints on LDM will be given 
for CTe, which plays the analogous role to the DM-nucleon scattering cross section, in (WIMP) 
DM scattering with nuclei. 

With these definitions, the cross section for a DM particle to excite an electron from level 1 to 
level 2 can be written as (see Appendix A.l) 

^ J ^s(^AEi^2 + ^ - q-v'^ X \Fj)Miq)\^\fi^2{q)f , (3.9) 

where fi^ 2 {q) is the atomic form factor for the excitation. It is given by 

fi^2iq) = j d^xil)*2{x)fi{x)e'-^'^ , (3.10) 

where 'fi and -02 are the normalized wavefunctions of the initial and final elecfron levels. We now 
apply fhis general resulf fo fhe special case of elecfrons in a periodic crysfal laffice, such as a semi¬ 
conductor. 


3.2.2 Excitation rate in a semiconductor crystal 

The periodic lattice of a semiconductor crystal has a continuum of electron energy levels, forming 
a complicated band structure (see Fig. 4). A small energy gap separates the occupied valence bands 
from the unoccupied conduction bands; exciting electrons across this bandgap creates mobile electron- 
hole pairs, which can be manipulated and detected. In order to perform practical calculations for this 
system, the true multi-body electron wavefunction must be replaced with a product of single-particle 
wavefunctions (this is a well-understood procedure, which we discuss further in §4). Once found, 
these single-particle wavefunctions can be used in Eqs. (3.9) and (3.10), giving the cross-section 
to excite an electron between specific energy levels. To find the total rate, these cross sections are 
integrated over initial and final electron levels, and over the DM velocity distribution. 


DM halo dependence. Neither the electron band structure, nor the electron wavefunctions, nor 
the DM velocity distribution are spherically symmetric. As noted in [9], the excitation rate will 
therefore depend on the orientation of the crystal with respect to the galaxy, an effect which may be 
extremely useful in verifying the DM origin of a signal. Here, however, we sidestep this complication 
by approximating the DM velocity distribution as being a spherically symmetric function gx{v)- All 
the relevant information about the DM velocity profile can then be encoded in the function r]{vrma) 
(see Appendix A.2), defined as 




dfv 

V 


Qxiv) Q{v - Vynin) 


(3.11) 


where 0 is the Heaviside step function. 
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When calculating rates, we assume a Maxwell-Boltzmann distribution with a sharp cutoff (we 
describe this in more detail, and give analytic formulas for ?y(nniin)> in Appendix B). The requirement 
of energy conservation is captured by r;min(Q, the minimum speed a DM particle requires in order 
for the electron to gain an energy with momentum transfer q (note that was also denoted as 
AEe in §3.1). This is given by 

Vminiq, Ee) = — + —^ . (3.12) 

q 2m^ 


DOS DOS 



Figure 4. Scissor corrected band structure for silicon (left) and germanium (right) as calculated with Quantum 
ESPRESSO [69] with a very fine k-point mesh. The horizontal dashed line indicates the top of the highest valence band. The 
four bands below the horizontal dashed line are the valence bands while the bands above the dashed line are the conduction 
bands. We also show the density-of-states (DOS) as a function of the energy for a very fine k-point mesh (blue) and for our 
243 k-point mesh (red). A Gaussian smearing of 0.15 eV was used to generate a smooth function. 


Differential rate. As we show in Appendix A.4, the differential electron scattering rate in a semi¬ 
conductor target (with the approximation of a spherically symmetric DM velocity distribution) can be 
written as 


(ii?crystal Py ,.7 _ 

,1 V, =— Acell 0-e a 

a in he 


X ^ / dlnq i —r]{vmm{q,Ee)) ) i^DM(g)^| /crystal (o',-Eg) 




(3.13) 


where py^ ~ 0.4 GeV/cm^ is the local DM density, E^ is the total energy deposited, and Aceii = 
Aftarget/Afceii is the number of unit cells in the crystal target. (Mceii = 2 x mce = 145.28 amu = 
135.33 GeV for germanium, and Mceii = 2 x msi = 56.18 amu = 52.33 GeV for silicon.) 
We have written this in such a way that the first line gives a rough estimate of the rate, about 
29 (11) events/kg/day for silicon (germanium) for py = 0.4 GeV/cm^, niy = 100 MeV, and cje — 
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3.6 X 10“^^ cm^ (the current limit from XENONIO [31]), while every factor in the second line is a 
roughly 0(1) number for the preferred values of q and E^. 

All the necessary details of the target’s electronic structure are contained in the dimensionless 
crystal form factor, /crystai(<?) Ef), which is a property purely of the target material and is independent 
of any DM physics. The computation of this form factor is one of the main results of this paper. 

Crystal form factor. In the periodic lattice of a semiconductor crystal, each electron energy level is 
labelled by a continuous wavevector k in the first Brillouin Zone (BZ), and by a discrete band index 
i. The wavefunctions of these states can be written in Bloch form, 

^ ^ u,(fc + , (3.14) 

G 

where the G’s are the reciprocal lattice vectors. Here V is the volume of the crystal, and the wave- 
functions are taken to be unit-normalized, so that 

Y,\ui{k + G)f = 1- (3.15) 

G 

Using this form for the wavefunctions, we can define the form factor for excitation from valence 
level {i k] to conduction level {i' k'}, 

f[ik,i'k',G'] ~ + G + G')ui{k + G). (3.16) 

G 

The crystal form factor required in Eq. (3.13) is then given by 

\f ( p m2 /■ Vce\\d?k Vcewd^k' 

|/crystal(g,f^e)| - ^ ^2vr)3 "" 

E, S(E, - + E.j) Y,Qi(<l-\i'-k + 0'|) ■ (3.17) 

S' 

(See Appendix A.4 for the derivation.) The band index i is summed over the filled energy bands, while 
i' is summed over unfilled bands, and the momentum integrals are over the 1st BZ. is the energy 
of level {ik}, and Uceii is the volume of the unit cell. The numerator in the first factor has units of 
energy, with value 27r^(amgUceii)~^ = 1-8 eV for germanium and 2.0 eV for silicon. The crystal form 
factor can be computed numerically using established solid-state computational techniques. Once it 
is known, it can be used to find event rates for any DM model and halo profile, using Eq. (3.13), along 
with Eqs. (3.7), (3.8), (3.11), and (3.12). We now turn to our own numerical evaluation of the crystal 
form factor. 

4 Numerical computation of the form factor 

Our aim is to compute the crystal form factor, given by Eq. (3.17), for silicon and germanium targets 
with low energy thresholds (;^ 30 eV). Once these are found, it is possible to calculate scattering rates 
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for any DM model. Calculating the form factor requires knowledge of the electron wavefunction co¬ 
efficients Ui{k + G) for all energetically accessible electron levels. To calculate these coefficients, we 
utilize the “plane wave self-consistent field” (PWscf) code within the Quantum ESPRESSO [69] 
package, based on the formalism of DFT. We then input these into our own postprocessing code, 
QEdark, to calculate the form factors. In this section, we summarize the key conceptual and numeri¬ 
cal details of our computation. We provide a review of DFT in Appendix F, detail the approximations 
used in the computation of the wavefunctions, and lay out the numerical methods. In Appendix C, we 
discuss the convergence of our computation. 

4.1 Computational framework 

It is impossible in practice to obtain the exact many-electron wavefunctions that describe interacting 
electrons in a many-body system such as a crystal. Instead, several methods exist to obtain excellent 
numerical approximations to these wavefunctions. We use DFT, which reformulates the interacting 
quantum many-body problem in terms of functionals of the particle density n(r). For the case of 
electrons, the Hohenberg-Kohn theorems [91] imply that all properties of the interacting system are 
determined once the ground-state electron density is known. In order to obtain the ground-state den¬ 
sity, we use the Kohn-Sham method [92] to map the system of interacting electrons into a system 
of independent electrons under the presence of an auxiliary potential that produces the same ground- 
state density. After this mapping, one has to solve the much simpler non-interacting electron system 
in order to obtain the ground-state energy and electron density. 

The mapping from an interacting to a non-interacting many-body system comes at the expense of 
having to use an approximate auxiliary potential. Typically this potential is split into the mean-field 
Hartree potential and an exchange-correlation potential. The latter captures the quantum mechani¬ 
cal effect of having identical electrons and also attempts to capture the correlation energy among the 
interacting electrons. The exchange-correlation potential is not known exactly and needs to be ap¬ 
proximated. We use the Perdew-Burke-Ernzerhof (PBE) functional [93], which belongs to the class 
of the Generalized Gradient Approximations (GGA). We discuss this further in Appendix E. 

Both silicon and germanium have a diamond lattice structure that contains two atoms in the unit 
cell. There are two s-shell and two p-shell valence electrons per atom (3s and 3p (2s and 2p) for 
germanium (silicon)), which makes a total of 8 electrons per cell. This translates to 4 valence bands, 
since each band is doubly degenerate in electron spin. In silicon, the core electrons have binding 
energies of ~100eV and above, and so are irrelevant for the energies we consider here. One must 
take more care with germanium, since the 3d electrons have binding energies of ~30eV, and so can 
be relevant for the higher energy thresholds we consider here.^ In the computation, energetically- 
inaccessible core electrons can be replaced with a pseudopotential, which increases the computational 
efficiency by reducing the number of initial states required, and by reducing the resolution needed to 
describe the wavefunctions. We use ultrasoft pseudopotentials [94] in place of all but the outer two 
s-shell and two p-shell valence electrons. Eor germanium, we also use a pseudopotential that allows 
us to treat the 3d electrons as valence states. As a result, the computational cost for germanium is 

"'We thank the authors of [68] for discussions regarding this point. 
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slightly higher than that of silicon. We use an empirical “scissor correction” approach [95, 96] to set 
the band gap to 1.11 eV for silicon and 0.67 eV for germanium [97].^ 


4.2 Discretization procedure and cutoff choices 

In order to obtain the crystal form factor with a hnite computation, several modihcations must be 
made to Eq. (3.17): 


• Binning in q and Eg. The form factor must be evaluated for hnite grid of q- and i7e-values. 
We do this by averaging over bins of equal width in q and Ef>'. 


I/. 


(binned) 

crystal 


(^n; Em 


)|2= 

~ Jqr^-lAq Aq 




.-|A£ 


dE' 


/crystal(9',^')|'- (4-1) 


Here q^ is the central value of the n*^ q bin, and E^ is the central value of the m**^ energy bin, 
and Aq and AE are the widths of the bins. We use 500 £'e-bins with AE = 0.1 eV and 900 
q-bins with Aq = 0.02 arrie- 


• Discretization in k. The continuum of fc-values in each energy band must be replaced with 
a discrete mesh of representative A:-points. The fc-integrals in Eq. (3.17) are then replaced with 
hnite sums: 



Here Vbz is the volume of the Brillouin Zone, V];eii is the volume of the crystal’s unit cell, and 
are the weightings of the k-points, with = 2 (following the convention of Quantum 

ESPRESSO). We use a uniform 243 /c -point mesh. 


• Cutoff in G, G. The wavefunctions are expanded in a hnite size plane-wave basis whose 
reciprocal lattice vectors satisfy the “kinetic energy” cutoff (really a cutoff in the space of G- 
vectors) 


2me 


< T'cut • 


(4.3) 


Note that since q= \k' — k + G'\, and since |Gmax| 1^1 and \k'\, the momentum transfer q 
essentially has a cutoff of \/2me77cut- We choose a value of Ecut = 70 Ry, which allows us to 
sample a large enough q space to obtain 0(1%) accuracy for our rate calculations. 


• Energy bands As discussed above, we consider initial electron states in the 4 valence bands 
for silicon and the 4 valence bands -i- 10 outer core bands (corresponding to the 3d-shell elec¬ 
trons) for germanium. We include hnal-state energy bands up to the 52°^ conduction band in 
both germanium and silicon. The lowest conduction states not included are about 57 eV above 
the band gap, while the highest energy core states not included are more than 60 eV below the 
band gap. Our choice of bands therefore fully covers any energy transition below ~57 eV. 

^These values are measured at 300 K, and change by 5-10% as the temperature approaches 0 K [94]. The effect of this 
on our results is a few percent and therefore negligible. 
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We can now write the form factor in the form that is implemented in our numerical code: 


I ^(numerical) / ^ n.| 2 _ 27r^ (amgl4ell) ^ ^k^k'\f |2 

Kcrystal \Qn^ m)\ ^ / j/ j / j 2 2 ,G']\ ^ 


0 1 - 


^i'k' 


i/c 


iAi? 


0 1 - 


\\k' -k + G'\-qn\ 


(4.4) 


Note that the first line here represents summing over bands, /c-points, and reciprocal lattice vectors, 
and calculating the contribution to the form factor from each. The sums are all over finite ranges as 
discussed above. The second line represents adding each contribution to the appropriate {q, Ef.} bin. 
We present the results of our computation, including prospects for upcoming experiments, in §6. In 
Appendix C we discuss convergence with respect to the choice of /c-point mesh and Ecut- 


5 Conversion from energy to ionization 

The calculation described in the previous two sections gives the DM-electron scattering rate in a 
semiconductor crystals as a function of the total energy deposited by the dark matter. Eg. 

However, experiments will not directly measure the deposited energy itself, but rather the ion¬ 
ization signal Q - i.e., the number of electron-hole pairs produced in an event. Linking the two is a 
complicated chain of secondary scattering processes, which rapidly redistribute the energy deposited 
in the initial scattering. 

A realistic treatment of the conversion from energy to ionization is a crucial step in calculating 
the sensitivity of experiments. Unfortunately, exact modeling of the secondary scattering processes is 
extremely challenging and is beyond the scope of this paper. Instead, we assume a linear response, 
which we believe does a reasonable job of capturing the true behavior. Specifically we assume fhaf, 
in addition fo the primary electron-hole pair produced by the initial scattering, one extra electron-hole 
pair is produced for every extra e of energy deposited above the band-gap energy. Here e is the mean 
energy per electron-hole pair as measured in high-energy recoils. The ionization Q is then given by 

Q{Ee) = 1 -|- [{Ee — Egi^p)/s\ , (5.1) 

where [xj rounds x down to the nearest integer, e and the band-gap energy ii'gap are measured to 
be [97, 98] 


f 3.6 eV (silicon) f 1.11 eV (silicon) 

£■ = N , -^gap = \ ■ (5-2) 

[2.9 eV (germanium) [0.67eV (germanium) 

We devote §5.1 and Appendix E to a discussion motivating this simple treatment. We emphasize 
that, while our treatment is approximate, it (a) is quite separate from the systematic, first-principles 
calculation of dR/dEe described in §3 and §4, and does not affect that calculation’s accuracy; (b) is 
probably conservative, since it does not account for fluctuations that could push a low-energy event 
above the ionization threshold; and (c) should be possible to improve upon in the future, both with 
better theoretical modeling and with experimental calibration. 
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5.1 Understanding the secondary scattering processes 


It is experimentally well-established that for high energy electron recoils keV), the ionization 
signal is directly proportional to the deposited energy, with a constant average energy e deposited per 
electron-hole pair created, 


(Q) == f ■ 


(53) 


e is several times the bandgap energy, accounting for the fact that only a fraction of the energy de¬ 
posited goes directly into pair production. Fluctuations around the average ionization are quite small, 
with the Fano factor (defined as the ratio of the variance to the mean) measured to be [99, 100] 


i7= 


(Q) 


0.1-0.15. 


(5.4) 


At the low energies we are interested in, 0(1-50 eV), the energy-ionization relationship has 
not been directly measured. Fortunately, there is reason to expect that the high-energy response can 
be extrapolated to lower energies. It has long been understood (see e.g. [98, 101]) that following a 
high-energy electron recoil, an electronic cascade occurs that rapidly redistributes the energy between 
many low-energy electrons and holes. Roughly speaking, any electron or hole is expected to re¬ 
scatter and create an additional electron-hole pair, so long as it has sufficient energy to do so. This 
repeats, distributing the energy over an exponentially increasing number of electron-hole pairs, until 
all electrons and holes have energy below the pair-creation threshold. Note that this threshold is larger 
than the band gap energy due to the constraints of momentum conservation [98]. The excess energy 
carried by the electrons and holes after the cascade is slowly lost to phonons, as is a fraction of the 
energy during the cascade. As a result of the cascade, the vast majority of secondary scatterings that 
occur after the initial electron recoil are low energy scatterings. This means that, for example, a single 
lOkeV electron recoil is approximately equivalent to 100 recoils with 100 eV each, or 1000 recoils 
with 10 eV each. This justifies fhe exfrapolafion of fhe high-energy behavior fo low energies. 

The linear response described by Eq. (5.1) is nof fhe only fracfable approach. Ofher, less sim- 
plisfic approaches can be faken wifhouf resorting fo a full firsf-principles freafmenf. For comparison, 
in Appendix E we consfrucf a phenomenological Monfe Carlo model of fhe secondary scaffering cas¬ 
cade, following [101]. The model is infended fo capfure fhe general fealures of fhe cascade, wifhouf 
knowledge of fhe specific microscopic strucfure of fhe fargef maferial. The model reproduces fhe 
known high-energy behavior well wifh only a single funable paramefer, and can be used insfead of 
Eq. (5.1) when calculating DM scaffering rafes. Unlike fhe linear freafmenf, fhe Monfe Carlo model 
predicfs flucfuafions abouf fhe mean, which can have an imporfanf effecl for DM masses fhaf are righf 
on fhe edge of defecfabilify. Eor typical masses, however, we find fhaf fhe fwo approaches agree fo 
wifhin a few lO’s of percenf (see Eig. 17). We conclude fhaf fhe linear frealmenf of Eq. (5.1) is a 
reasonably realisfic approximafion, and posfpone a more careful frealmenf fo fulure work. 


6 Results 

In Ihis secfion, we presenl fhe resulls of our calculalion of fhe DM-eleclron scattering rales in sili¬ 
con and germanium defeclors. We show fhe polenlial reach for single-eleclron-sensilive experimenls. 
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Figure 5. The crystal form factor |/crystai(g, Ee)'\^ as a function of q and for silicon (left) and germanium (right) 
(see Eq. (3.17)). In the region below the solid line, Umin > Uesc + ve for any DM mass, and electron scattering is thus 
kinematically inaccessible. The dashed line corresponds to Umin = 300 km/s (a typical DM halo velocity) in the heavy 
DM limit; the region below this line is only kinematically accessible to DM particles with velocities larger than the average 
velocity. For energies above ~ 10 eV, the scattering rate is suppressed by both the form factor and DM velocity distribution. 
We see that the 3d electrons in germanium give a sizable contribution to |/crystal(q, Ee)\‘^ for Ee > 25 eV. 


as well as the effect of higher experimental thresholds. We also give the full recoil spectra and the 
annual modulation fraction, which may be crucial for discriminating a possible signal from back¬ 
ground. Lastly we discuss near-term prospects, focussing on upcoming searches expected from the 
SuperCDMS and DAMIC collaborations. 

Experimental thresholds are set in terms of the ionization signal Q (the number of electron-hole 
pairs produced in an event) rather than the deposited energy Ee- In the following results, we model 
the conversion of deposited energy to ionization with the linear treatment described in §5. We take 
the DM halo to have a local density of pdm = 0.4 GeV/cm^ [102, 103], and a Maxwell-Boltzmann 
velocity distribution with a mean velocity vq = 230 km/s and escape velocity Vesc = 600 km/s, and 
we take the average Earth velocity to be r;E = 240 km/s (see Appendix B for explicit formulae). In 
Appendix C we discuss the numerical convergence of our results. 

Event rates as a function of Q, for an extensive range of DM masses, are available online at this 
link. The crystal form-factor, as a function of q and E^, is also available there. Using Eq. (3.13), the 
information online can be used to re-derive rates using a different DM form-factor or velocity profile, 
or using a different treatment of the energy-to-ionization conversion. 

6.1 The crystal form factor 

Much of the behavior of the scattering rates can be understood from the behavior of the crystal form 
factor, I /crystal (<?5 in Eq- (3.17). We show the crystal form factor in Eig. 5 as a function of q 


- 19 - 














[MeV] 



Figure 6. Dark Matter-Electron Cross-Section Sensi¬ 
tivity: The 95% C.L. exclusion reach in the DM-electron 
scattering cross-section, of an experiment with 1 kg- 
year exposure and zero background events, for different ex¬ 
perimental thresholds. Solid (dashed) lines show the reach 
for silicon (germanium) targets. Ionization thresholds of 1, 
5, and 10 electron-hole pairs are shown with blue, green, 
and red lines, respectively. The corresponding energy 
thresholds are 0, 11.6, and 26.1 eV in germanium, and 0, 
14.4, and 32.4 eV in silicon. The gray shaded region shows 
the existing constraint fromXENONlO data [31]. The three 
plots assume different DM form factors, Fumiq) = 1, 
arrielq, (am^jqf', corresponding to different DM mod¬ 
els. 


and £'e, for both silicon and germanium. The rapid fall-off as q increases is clearly visible. The solid 
line in the figure corresponds to V rr,\n = v^sc + ve from Eq. (3.12) as —)■ oo. The region below 
this line is kinematically inaccessible for any DM mass. The dashed line uses the velocity of a typical 
DM particle in the halo, i.e. fmin = 300 km/s. We see that larger recoil energies require larger q, for 
which the crystal form factor is suppressed. The implication of this is that the DM-electron scattering 
rates increase dramatically for smaller recoil energies, resulting in a dramatic increase in sensitivity 
as detector thresholds are lowered. 
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Figure 7. Spectrum of events as function of the ionization signal Q. The different lines show different DM masses and 
form-factors, as indicated, in a silicon (left) or germanium (right) target. The spectra are normalized to 1 in the first hin. 
The top axes show the values of the deposited energy Ee corresponding to the edges of the bins. Since the typical Ee is 
around a few eV, the distributions peak at Q = 2 (see §3.1). 


6.2 Cross-section reach versus detector threshold 

In Fig. 6, we show the sensitivity to the DM-electron scattering cross section, ae, versus the DM 
mass, m-^, for hypothetical silicon- and germanium-based experiments with a 1 kg-year exposure 
and zero background, and with various detector thresholds. The curves show 95% C.L. limits, i.e. 
3.6 signal events. The blue, green, and red lines show ionization thresholds, Qth, of 1, 5, and 10 
detected electron-hole pairs, respectively, which correspond to deposited energies, E^, of 0.67, 12.3, 
and 26.8 eV in germanium, and 1.1, 15.5, and 33.5 eV in silicon (to get the corresponding ionization 
energy thresholds, subtract 0.67 eV for germanium and 1.1 eV for silicon from these numbers, see 
Eq. (5.1)). The three plots show results for different DM form factors, corresponding to different 
classes of DM models: Fdm(<?) = 1 {top left), Fj}m{q) = aruelq (top right), and Fj^uig) = 
{amelq)"^ {bottom), see §2 for details. 

As expected, the reach dramatically improves when the threshold is lowered, since the crystal 
form factor strongly suppresses the rate for high electron recoil energies. This improvement is most 
pronounced for Fdm(9) = {ame/q)"^, since lower q tends to correspond to lower recoil energies. 
With a single-electron threshold, the difference in sensitivity for silicon and germanium targets can 
be accounted for by the fact that germanium is 2.6 times heavier, and so has correspondingly fewer 
valence electrons per kg. However, germanium targets are sensitive to slightly lower DM masses 
due to their lower band-gap. In addition, germanium targets become comparably more sensitive than 
silicon targets for ionization thresholds of Qth 9 due to the additional contribution from the 3d-shell 
electrons (see below). 
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Fig. 7 shows the spectrum of events as a function of the ionization signal Q, for different DM 
form-factors and masses, in silicon (left) and germanium {right) targets. The fast fall-off with Q shows 
the large gain to be made from lowering experimental thresholds towards a 1- or 2-electron threshold, 
especially for the steeper DM form-factors and for lower DM masses. The shape of the spectra may 
be useful in discriminating a signal from background. 

In germanium, the 3d-shell electrons dominate the rate for > 24 eV, corresponding to Q > 9 
electron-hole pairs (the 3d-shell electrons lie about ~15 eV below the bottom of the valence band and 
~24 eV below the bottom of the conduction band). The intuitive reason for the 3d-shell electrons 
dominating over the valence electrons at large can be seen from Eqs. (3.5) and (3.6). The typical 
velocity of the 3d-shell electrons, and hence the typical momentum transferred from the DM, is larger 
than for the valence electrons, so that the 3d-shell electrons can dominate if they are kinematically 
accessible.® We show the effect of neglecting the 3d-shell electrons in Fig. 15 in Appendix D, where 
we compare the cross-section reach and the recoil spectrum generated by DM scattering with and 
without the inclusion of the 3d-shell electrons. The effect is significant for ionization thresholds 
above Qth ~ 7 or 8, but not important at lower thresholds. 

We note that there are some differences between our results and those in [9, 25]. In [9], only the 
case Qth = 1 was considered and we find fhaf fhe new compufafion predicfs a somewhaf lower rale. 
We find lhal fhe shape of fhe recoil speclra in [25] is noticeably differenl from ours, which gives rise fo 
several differences in fhe expecled limils for fhe differenl Q Ihresholds. Furlhermore, for germanium, 
we also include fhe 3d-shell elections, which can be imporlanl as discussed above. 

6.3 Comparison with existing XENONIO limit and discussion of background 

We see from Fig. 6 that to surpass the existing limits obtained with XENONIO data [31], agermanium- 
or silicon-based experiment with an ionization threshold of 10 elections would require a background- 
free exposure of less than 1 kg-year. However, with a single-election threshold, such an experiment 
would surpass the XENONIO limit at all masses with a background-free exposure of around 1 kg-day 
for Fdm( 9) = 1, lOg-days for FbM(9) = aruelq, or just a 1 g-day for Fdm(9) = {anielq)'^. In 
addition, with any exposure such an experiment would place the first ever bounds in the ~l-5MeV 
mass range, below the threshold of the XENONIO search. The XENONIO detector had a threshold 
of one electron with an 0(1) detection efficiency, but to obtain one electron required an energy of 
at least 12.4 eV to overcome the binding energy of an electron in the outer shell. ^ Moreover, the 
background in the XENONIO data was much larger than conventional nuclear-recoil background, so 
that the number of DM events leading to single (two, and three) elections was only limited, at 90% 
C.E., to be less than 8,550 (1,550, and 330) counts/kg/year, respectively. 

While the single-electron background in the XENONIO data was rather large, its origin is likely 
specific to its dual-phase detector setup. Many of the single electron backgrounds likely had one, 
or a combination, of the following origins [31]: (i) electrons, trapped in the potential barrier at the 

*For even larger deposited energies, Ee > 100 eV, it is likely that the deeper shells will dominate the rate, for both 
silicon and germanium. 

’Note that in [31], the electrons were treated as bounded inside free (xenon) atoms, unlike the electrons in the semicon¬ 
ductor targets here. 
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liquid-gas interface, were randomly drawn into the gas phase (these transiently trapped electrons likely 
originated from other background events that caused xenon atoms to ionize); (ii) photo-dissociation of 
a negatively charged 02 -ion, which received its negative charge from a drifting electron that originated 
from another event; (iii) field emission from the cathode. XENONIOO and LUX may face similar 
challenges, although an analysis is still in progress. 

The semi-conductor targets will not suffer from these same detector-specific backgrounds. They 
will, of course, have fheir own unique experimenfal challenges fo deal wifh, including defector noise 
and dark currenf, as we will discuss in more defail in §6.5 for DAMIC and SuperCDMS. These will 
likely be fhe limiting insfrumenfal facfors in selling fhe fhreshold for a parlicular experimenf. Once 
Ihese challenges are overcome, one needs fo deal wifh fhe physics backgrounds. As argued in [9], 
neulrinos are nof an imporlanl source of background even for fhe largesl exposures considered in fhis 
paper (0(20 kg-years) for SuperCDMS, see §6.5). Compton scalfering or olher evenls lhal produce 
recoiling eleclrons will usually lead fo a much larger energy deposition and mosl of Ihem could fhus 
be vetoed, allhough some backgrounds will persisl fo fhe lowesl energies. The size of fhis background 
will depend on fhe shielding and purify of fhe materials around fhe deleclor; for SuperCDMS al SNO- 
LAB, fhe Compton background is estimated to be 0(6 x 10“^) events/kg/day/keV [104]. Assuming 
that it is flat at low energies, this translates into 0(0.04) events/eV for 20 kg-years, which would be 
negligible. We do not expect there to be backgrounds from neutrons, and (cosmogenic) x-ray tines 
will tie well above our energies of interest. Surface events and other, unknown, backgrounds may be 
present at low energies. As experiments reach the required sensitivity to probe the few-electron events 
expected from LDM scattering off electrons, a better understanding of all backgrounds will emerge 
allowing for an attempt to mitigate them if necessary. It is possible that a spectral analysis of a signal 
will further allow for the removal of some background events. Our assumption of zero background 
for the plots should be taken as the best-case scenario. 

6.4 Annual Modulation 

Even with a significant background event rate, it may be possible to distinguish a signal from back¬ 
ground using annual modulation, as long as the background is stable on year time-scales. Annual 
modulation is a distinguishing feature of a LDM scattering signal [9, 105], occurring due to the change 
in the earth’s velocity through the DM halo as it rotates around the Sun. Eor a standard smooth and 
isotropic DM velocity distribution, the modulation is approximately sinusoidal with year period and 
a peak around June 2nd (the presence of DM streams or non-trivial DM structure may complicate 
this, as may gravitational focusing by the Sun [68, 106], which we do not include). The modulation 
fraction, /mod^ is defined fo be fhe ratio of fhe amplifude of fhe modulating signal fo fhe median signal 
rate, which for our assumed halo profile (see Appendix B) corresponds fo 

/n.„d = R„ = 2 = flSept 2 ( 6 - 1 ) 

ZKq 

where i?j represenfs fhe rate af fime of year i. 

Eor DM scaffering off elections, fhe modulafion fraction can be significanlly larger fhan for fhe 
usual elastic scaffering of (heavy) WIMPs off nuclei. As we saw in Eig. 5 (see discussion in §3.1), 
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Figure 8. Annual modulation fraction of the DM signal, /mod. Left: /mod as a function of the ionization threshold 
Qth, in silicon, for 10 MeV (blue) and 1 GeV (black) DM masses. The solid (dashed) lines correspond to Fum = 1 
(_Fdm = {anie/q)^)- The top axis indicates the energies corresponding to the edges of the bins. Center: Same as left 
for germanium. Right: /mod as a function of DM mass, for ionization thresholds of 1 and 5 electron-hole pairs. The 
solid (dashed) lines correspond to Tdm = 1 (Tdm = {anielq)^), while the blue (green) lines correspond to germanium 
(silicon). 


DM-electron scattering relies on the tail of the DM velocity distribution, especially for energies above 
~ 5 -10 eV. We plot the modulation fraction in Fig. 8. The left and center plots show /mod as a function 
of ionization Q for different masses and DM form-factors for the two elements, /mod rises from a few 
percent for single-electron events to above 10% for events with more than ~ 3 electrons. Comparing 
with the spectrum in Fig. 7, we see that there is a trade-off between modulation fraction and event rate. 
Events with several electron-hole pairs provide large modulation without sacrificing too much in the 
rate, and may give the best prospects for annual modulation searches depending on the background. 
The modulation fraction also rises near the mass threshold, as we show on the right of Fig. 8 for 
ionization thresholds of (5th = 1 and 5. Note that the high-mass value of /mod for tho single-electron 
threshold, at 4 -6%, is larger than the values in the Q = I bin of the left and center plots, because the 
total rate is not dominated by the single-electron events. 

Once a signal is found in an electron scattering search, increasing the exposure of the experiment 
until the annual modulation can be tested will be a crucial step in claiming a DM discovery. In Fig. 9 
we show the her discovery reach of an annual modulation search in the mass-cross-section plane. We 
calculate this cross section by requiring 

^S/^Stot + B = 5 , (6.2) 

where AS = fmodStot is the modulation amplitude, Stot is the total number of signal events, and 
B is the number of background events. The thick curves in Fig. 9 show the discovery reach for 
different thresholds and DM form-factors, assuming a background-free exposure of 1 kg-year. (A 
non-zero background will of course weaken the reach, following the equation above.) This figure 
mirrors Fig. 6, which shows the exclusion reach obtained using a simple counting search instead of 
a modulation search, but otherwise with the same assumptions. The curves of Fig. 6 are replotted as 
the thin curves in Fig. 9, for comparison. We see that a substantial discovery reach is possible with 
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Figure 9. Annual Modulation: Discovery Reach: The 

5a discovery reach in the m^-ae plane of an annual mod¬ 
ulation search for DM-electron scattering for an experi¬ 
ment with a 1 kg-year background-free exposure. Solid 
(dashed) lines show the reach for silicon (germanium) tar¬ 
gets. Ionization thresholds of 1, 5, and 10 electron-hole 
pairs are shown with blue, green, and red lines, respec¬ 
tively. The corresponding energy thresholds are 0, 12.3, and 
26.8 eV in germanium, and 0, 15.5, and 33.5 eV in silicon. 
The gray shaded region shows the existing constraint from 
XENONIO data [31]. The three plots assume different DM 
form factors, as indicated, corresponding to different DM 
models. Thin lines are from Fig. 6, showing the exclusion 
reach of a search with the same exposure seeing no events. 


a 1 kg-year exposure: at low masses for Fdm(^) = and at all masses for Fdm(9) = arnejq or 

{anielq)^. 

Finally, we comment that taking into account the directional (sub-daily) modulation, which is 
expected in crystalline detectors, will further allow for an improved sensitivity to a DM signal. As 
discussed in §3.2.2, we have averaged-out such directional effects in this work, and we postpone their 
study to future work. 

6.5 Prospects for near-term experiment 

In this subsection, we discuss the near-term prospects for electron-scattering searches with the DAMIC 
and SuperCDMS experiments. 
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6.5.1 DAMIC 


DAMIC [107-109] uses thick, fully-depleted silicon CCDs for their target material. These CCDs are 
ten times more massive than conventional CCDs, allowing them to be competitive targets for DM 
direct detection. In [108], DAMIC used one 0.5 g CCD to perform an engineering run, obtaining an 
exposure of 107 g-days. They were able to constrain DM-nncZear scattering for DM masses almost 
as low as 1 GeV. Work is ongoing to increase the total mass of the detector (by using more CCDs) as 
well as the detector’s sensitivity to low threshold energies (by using so-called “Skipper CCDs”) [67]. 

The first direct detection limit using a semi-conductor target. Here we investigate the (albeit 
weak) constraints on DM-electron scattering from their existing result, and give reach estimates based 
on their projected detector improvements. For the engineering run [108], DAMIC used a single 0.5 g 
CCD, for an exposure of 107 g-days. They obtained the following values for the read-out noise and 
the dark current: 

(i) A readout noise of below 2 electrons/pixel, corresponding to 2 x 3.6 eV = 7.2 eV of r.m.s. read¬ 
out noise. The CCD has about 4.2 million pixels, so that one requires a threshold of ~ 13 
electrons (~ 47 eV) for the noise to produce a signal above threshold in less than one pixel. 
DAMIC chose a threshold of ~ 40 eV (~ 11 electrons) for the search for DM-nuclear scatter¬ 
ing; we expect ~ 35 pixels to reach this threshold. In our recast of their data for DM-electron 
scattering below, we will use the same 40 eV threshold. 

(ii) A dark current of ~ 1 electron/CCD/day (at the chosen 120 K operating temperature). Since the 
exposure of the CCD is a few hours, before being read-out within a few minutes, the threshold 
is limited by the read-out noise, and not the dark current. 

We can use the result in [108] to constrain DM-electron scattering. We will assume that the efficiency 
to select electron recoil events is the same as selecting nuclear recoils, i.e. 7 x 10“^. Fig. 12 in [108] 
shows the data that was recorded by DAMIC, and we see that zero events were recorded in the first bin 
above the threshold (40 eV to 100 eV). This may be as a result of the efficiency being very low in this 
bin; nevertheless, it could also be a sign that backgrounds may be small at such low energies, boding 
well for future runs with even lower energy thresholds. In any case, this information is sufficient to 
derive the current DAMIC limit on LDM, which we show with a green shaded region (bounded by a 
green line) in Fig. 1 . We see that with the current threshold, the form-factor suppression is too large 
for this constraint to compete with the existing XENONlO-based limit [31]. Nevertheless, this is the 
first direct detection limit for sub-GeV DM using a semi-conductor target. 

Projections for future DAMIC runs with improved “Skipper” CCDs. There are two main chal¬ 
lenges that need to be overcome by DAMIC (and similar experiments) to allow them to push to low 
thresholds [67] : (i) reduce the noise in reading out the ionization deposited in the detector, and (ii) 
reduce the dark current. The read-out noise can be reduced substantially by taking more time to read 
the CCD, while the dark current (i.e. genuine electron-hole pairs produced by thermal excitations in 
the silicon substrate) can be reduced by lowering the temperature and improving the quality of the sil¬ 
icon. The contribution from the dark current will increase with the readout time, so it will take some 
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Ionization fhreshold 

dark currenf: 5x10 ^elec./pixel/day 

dark currenf: 10 ^elec./pixel/day 

1 kg-day 

100 g-year 

1 kg-day 

100 g-year 

Qth = 1 

1.6 X 10^ 

5.8 X 10® 

320.0 

1.2 X lO'^ 

2 

1.7 X 10^ 

6.1 X 10"^ 

6.7 X 10"'^ 

2.4 X 10'^ (*, t) 

3 

0.1 (A) 

4.2 (o) 

9.3 X 10"^® 

3.4 X 10"^^ 

4 

6.0 X 10"® 

2.2 X 10"^ 

9.6 X 10"25 

3.5 X 10"23 


Table 1. Expected number of events with at least Qth electron-hole pairs under different assumptions for the dark current 
and exposure (5 x 10“®electrons/pixel/day and 10“^electrons/pixel/day) and assuming either (i) 4 CCDs (10 g) and an 
exposure of 1 kg-day; or (ii) 40 CCDs (100 g) and an exposure of 100 g-years. In both cases, we assume that it takes 
one hour to read the entire CCD, and that it is read continuously. Projected exclusion reach based on a simple counting 
experiment are given in Fig. 1 for the entries marked with a single star (*). A projected discovery reach based on seeing the 
annual modulation of the signal, with negligible background, is also given in Fig. 1 for the entry marked with a dagger (f) 
(we find that the prospects from annual modulation for the entry marked with a diamond (o) are very similar). See text for 
details. 


optimization to find a way to reduce the readout noise while keeping the contribution from the dark 
current manageable. Lowering the temperature also reduces the electron mobility in the substrate, 
requiring a careful trade-off. Here we project what future data runs can achieve with the improved 
DAMIC Skipper CCDs. 

The DAMIC Collaboration has been working on so-called “Skipper CCDs”, which will reduce 
the r.m.s. read-out noise down to 0.2 electrons/pixel/day, with the possibility of going down to 0.1 elec- 
trons/pixel/day [67, 108]. This is done with a new output circuit that enables multiple read-outs. The 
size of the CCDs can be anything up to 4 x 4 = 16 million pixels (Mpix) [110], but we will assume 
8 Mpix for the projections below. The 0.2 (0.1) electrons/pixel/day correspond to read-out noise of 
0.72 (0.36) eV; assuming gaussian noise, the 0.1 electron/pixel/day will allow for sensitivity down 
to single electrons, although the 0.2 electrons/pixel/day may require a threshold of two electrons to 
avoid the noise faking a DM signal. However, for such low read-out noise, the dark current becomes 
the limiting factor in determining the energy threshold. Significant non-gaussian tails could change 
this conclusion. 

The dark current has been measured currently at 5 x 10“^ electrons/pixel/day [110]. As men¬ 
tioned above, lowering the operating temperature and improving the silicon substrate quality will de¬ 
crease this, and it is reasonable to expect further improvement; the theoretical lower limit is 0(10“^) [110]. 
Below we provide projections under both assumptions. 

The effect of the dark current on the threshold depends on the number of pixels and the exposure 
length of the CCD. The CCD is read pixel-by-pixel, and can be read continuously from one side to 
the other, before cycling back again to the beginning. We will assume that the 8 million pixels of 
the CCD are all read in one hour, so that its exposure is one hour for the purposes of calculating the 
dark current. We will consider the following two scenarios for the number of CCDs, the mass, and 
exposure (we assume an efficiency of 1 for making our projecfions below): 

(i) There are currenfly four profofype skipper CCDs, each wifh a mass of 2.5 g, which were pro- 
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duced as part of an R&D project (these will be deployed at the MINOS near site this year). For 
our first set of projections, we will assume that data is taken over 100 days (livetime), for a total 
exposure of Ikg-day. 

(ii) If the testing of the skipper CCDs at MINOS goes well, one can expect that several more of 
them will be deployed to search for DM. Thus, for our second set of projections, we will assume 
that 40 CCDs are deployed (for a total mass of 100 g) and that data is taken again over 365 days 
(livetime) for a total exposure of 100 g-years. 

Table 1 gives the expected number of events with at least Q electron-hole pairs for the two scenarios 
assuming Poisson statistics. We see that DAMIC could have a threshold of 2 electrons if the dark 
current can be reduced to below 10“^electrons/pixel/day, but a 3-electron threshold is required for 
the present dark current rate of 5 x 10“^electrons/pixel/day. In Fig. 1, we show solid green lines 
that indicate the 95% C.L. prospects for the entries marked with a star (x), i.e. we show the cross 
section to obtain 3.6 signal events, assuming zero background events. We also show the reach of an 
annual modulation search for the entry marked with a dagger (f), i.e. a 2-electron threshold with no 
background in a 100 g-year exposure. This is shown by the dashed green line in Fig. 1. We checked 
that the prospects from annual modulation for the entry marked with diamond (o) are very similar. We 
see from these projections that DAMIC can significantly improve upon the current XENONIO limit, 
especially at the lowest DM masses. 

6.5.2 (Super)CDMS 

The CDMS experiment uses cryogenic solid state detectors operated at temperatures below 100 mK. 
In the WIMP search, they distinguish electron from nuclear recoils by measuring the ratio of the 
ionization versus phonon energy deposited into the crystal. This ratio will be smaller for nuclear 
recoils than for electron recoils. Here, we are interested in their ability to detect electron recoils. 

The signal from a low-energy recoiling electron can be dramatically enhanced by applying a 
relatively large bias voltage, H ~ 0(50 — 100 V), across the target material. The work done in 
drifting an electron-hole pair out of the crystal, el4, is emitted as Luke-Neganov [111-113] phonons, 
which will be picked up by the phonon sensors. This was done for “CDMSlite”, which has yielded 
electron recoil thresholds with an 0(1) detection efficiency of 170eV (i.e. 0(50) electrons) [89]. 
Here we discuss the prospects of future versions of SuperCDMS. 

In Fig. 1, we show 3 projections for SuperCDMS, two for silicon (with an exposure ~10 kg- 
years) assuming an electron-hole pair threshold, (5th^ of either 4 or 1, and a signal detection efficiency 
of 0.7, and one for an annual modulation search assuming Qth = 2 with the same exposure and 
efficiency. The cross-section reach for germanium is very similar to those of silicon. 

The electron-hole-pair thresholds are based on the following assumptions. The Qth = 4 thresh¬ 
old is based on the numbers used by SuperCDMS for Snowmass [63], while the Qth = 1 threshold is 
based on an ambitious but achievable best-case scenario. For Snowmass, a phonon energy resolution 
of 50 eV was assumed. As there may be non-gaussian tails, a 7a threshold was assumed, correspond¬ 
ing to a threshold of E^, = 350 eV. Taking the bias voltage to be 100 V, this translates into Qth = 3.5, 
which we round up to 4. For the second set of projections, we assume that further R&D can push the 
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noise threshold down to better than ~14 eV, which is ambitious but achievable in principle [114]. A 
7a threshold corresponds to ~100 eV, so that a bias voltage of 100 V is sufficient to achieve sensitivity 
to Qth = 1- In practice, the bias voltage can be optimized as well. A larger bias voltage would allow 
for a reduced threshold in terms of the number of electron-hole pairs needed to pass the phonon energy 
threshold, but could lead to breakdowns of the substrate. However, it has been demonstrated that the 
bias voltages needed for sensitivity to Qth = 1 are achievable in both silicon and germanium [66]. 

As can be seen from Fig. 1, SuperCDMS has the potential to improve drastically upon the existing 
XENONIO limit, especially at low DM masses. 

7 Conclusions 

Direct detection experiments have so far primarily focused on searching for WIMPs, and as a result 
of an intense research effort, the path forward in this direction is rather well-defined. Wifhin fhe 
nexf decade, WIMPs will eifher be found or become significanfly less mofivafed. However, ofher 
fheorefically mofivafed candidafes exisf fhaf could consfifufe fhe DM in our Universe. In fhis work, 
we focused on a class of DM candidafes fhaf have a mass befween a few-hundred keV fo a GeV. We 
showed fhaf fremendous progress can be made in exploring fhe direcf-defecfion paramefer space of 
fhese candidafes over fhe nexf few years, by searching for DM-induced elecfron recoils in experimenfs 
wifh fargefs fhaf consisf of semiconducfor maferials. The fechnology currenfly used in WIMP searches 
can be adapfed for such lighf-DM searches by improving fhe ionization sensifivify, and fhis is being 
acfively pursued. The backgrounds are expecfed fo be quife differenf in nafure fo fhose in WIMP 
searches, and fhere is reason fo believe fhaf fhey will be small and confrollable. 

The calculation of fhe DM-elecfron scattering rafe and fhe subsequenf elecfron recoil specfrum 
in semiconducfor fargefs is much more challenging fhan for DM-nuclei scattering. We have provided 
defailed formulae for fhe scattering rafe and recoil specfrum, expressed in ferms of simple DM prop¬ 
erties and a fargef-dependenf “crysfal form factor”, which encodes fhe quanfum sfrucfure of fhe fargef 
elections. We numerically calculafed fhe crysfal form facfor for germanium and silicon wifh our code 
QEdark, which is based on fhe soffware package Quantum Espresso fhaf calculafes fhe crys¬ 
fal wave functions and energy levels. Convergence fesfs indicafe fhaf our resulfs are accurafe af fhe 
few percenf level. QEdark will be publicly available af ddldm. physics . sunysb . edu, fogefher 
wifh fhe crysfal form facfors. Wifh fhese, upcoming experimenfs can derive fheir own sensifivifies or 
limits. 

The crystal form factor is a steeply falling function of the electron recoil energy. Consequently, 
even a small improvement in an experiment’s detector threshold translates into a significant increase 
in the sensitivity to DM-electron scattering. We have provided the projected sensitivity for a variety 
of experimental thresholds, showing that upcoming experiments including DAMIC and SuperCDMS 
can probe orders of magnitude of unexplored DM parameter space in the next few years. In addition to 
setting limits, sub-GeV dark matter can be discovered via its expected modulation signal. We showed 
that in the case of electron-scattering, annual modulation is sizable and could provide the necessary 
signal for discovery. Additional sub-daily modulation is expected due to the orientation-dependent 
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nature of scattering in crystalline detectors. We have ignored directionality in this work, deferring it 
to future study. 

Calculating the experimentally observable signal requires a conversion from energy deposition to 
the ultimate ionization signal. This conversion requires a detailed knowledge of the secondary scat¬ 
tering processes in crystals, at energies below the existing experimental sensitivity. We therefore used 
a phenomenological model for secondary interactions, and studied its possible systematic uncertain¬ 
ties using a Monte Carlo model. We find that our predictions suffer from systematic uncertainties of 
order a few tens of percent, and is likely conservative. Further theoretical and experimental study of 
secondary interactions would be useful to improve the modeling of this conversion. 

To summarize, our work provides the necessary tools for experiments which use semiconductor 
targets to search for sub-GeV dark matter to derive accurate limits. Technologies adapted from WIMP 
searches and currently under development can be employed in searches for sub-GeV dark matter. This 
highly-motivated direction in dark matter searches is a natural progression from the WIMP program, 
and we expect that it will take a leading role in the search for dark matter. 

Note added: 

While this work was being completed. Ref. [68] appeared, which also deals with DM-electron scatter¬ 
ing in germanium. Ref. [68] is complementary to our work, its main point being the effect of “grav¬ 
itational focusing” on the modulation signal of DM-electron scattering. Ref. [? ] derives scattering 
rates using a semi-analytic approach, which builds on the method of Ref. [25], but is significantly less 
detailed than the method we have presented here. We find thaf fheir results are comparable to ours 
within a factor of a few, but with some notable differences. In particular. Ref. [68] finds increasingly 
higher rates than us at increasingly higher recoil energies. Most strikingly. Ref. [68] finds that scat¬ 
tering of the 3d shell electrons dominates the total rate when it is kinematically accessible (we find 
the the 3d shells cause a bump in the spectrum, but with a rate subdominant to lower energy events). 
We attribute these differences to the inherent sensitivity of the calculation to the tails of the electron 
wavefunctions for energies above 0(10 eV), as we discussed in §3.1. 
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A Derivation of scattering rate formulae 


A.l General formula for dark matter-induced electronic transitions 


If a DM parficle scaffers wifh an elecfron in a sfafionary bound slate such as in an alom, if can excite 
fhe elecfron from an initial energy level 1 lo an excifed energy level 2, by Iransferring lo if energy 
A£'i _^2 and momenlum q. The cross secfion for Ibis process can be derived in a slandard way using 
non-relalivislic quanlum mechanics, bul here we derive if slarling from fhe usual formula for fhe 
cross-seclion in field-lheory, in order lo make easier connecfion wifh fhe underlying parficle physics. 
We Ireal fhe elecfron as being bound in a sialic background polenlial - in olher words we approximate 
fhe aloms as being infinilely heavy objecls which can absorb momenlum wilhoul recoiling. This is 
an excellenl approximalion (< 1% error), since fhe momenlum-lransfers we will be interested are 
lypically of order keV. 

The cross secfion for free 2 —)■ 2 scallering is given by 




1 

4^ 


d^q d^k' 1 
(27r)3 (27r)3 4:E^Ee 


{2TTf5{Ei-Ef)5^{k + q- 


k')\Miree{q)? , 


(A.l) 


where Ad free is the usual field-theory matrix element, and |Af p represents its absolute square aver¬ 
aged over initial spins and summed over final spins. 

If fhe elecfron were unbound, fhe non-relalivislic scallering amplifude would be given by 


Ap|d7int|Xp',efc) = CMiree{q) X {27rf5^{k - q- k '), (A.2) 

where eg) is plane-wave slate for a DM parficle of momenlum p and an elecfron of momenlum 
k, iFint is fhe inleraclion Hamillonian, and C is an unimporlanl coefficienl. However, because fhe 
election is bound if is instead given by 


(Xp>_g-,e2|d^int|Xp>,ei) = 


( 2^^3 Mk ){Xpnej;,\ 




int 




(27r) 


/ 1/^/3 h ~ ^ ~ ^ 

+ q)'4’i{k) , 


(A.3) 


where il)i, y )2 are fhe (unil normalized) momenfum-space wavefunclions of fhe initial and final elec¬ 
tion levels. We have purposefully used plane-wave normalization for bolh fhe free and bound election 
slates: (eg|eg) = (ei|ei) = (27r)3(53(()) = V, where V is fhe volume of space (which always cancels 
in fhe end). 
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To find the cross section for this excitation process, we can use the free 2 —)• 2 scattering cross 
section formula but with two replacements: one to account for the modified scattering amplitude, and 
the other to account for the different final-state phase space. Squaring Eqs. (A.2, A.3), we see that the 
bound-state scattering amplitude is accounted for by making the replacement 

V{27:f5\k - q- k')\Mi,ee? lA^freeP X V^\h^ 2 {q)? , (A.4) 


where fi^ 2 {q) is the atomic form factor. 




dfk 


fl{k + q)i)i{k). 


(A.5) 


j (27r)^ 

Fourier transforming Eq. (A.5) gives the definition given in Eq. (3.10). 

Since there is only one final electron state being considered, we also need to remove the usual 
final-state phase space integral: 

^ (fk' 


free-electron phase space = V 


(2vr) 


1 


(A.6) 


Combining Eqs. (A.l, A.4, A.6), we can write the formula for the cross-section for a DM particle 
to excite an electron from level 1 to level 2: 


1 


0-Vl^2 = 


dfq 


1 


4.E' E'^ J (27r)3 4.E^E, 


-27:5{Ei-Ef)\Mi,ee{qW x \h^2{q)Y 


Since we are in the non-relativistic regime, the energies are given by 

1 2 

Ei = my. + me + -m^v + Ee^i 

„ \m^v — q P „ 

Ef = m^ + me + -—-h Ee,2 ■ 

2m^ 

Using the following definitions [9] to parametrize the underlying DM-electron coupling 

|Adfree(9)P = | At free (a^e) ^ X |Fdm(Q')|^ 




/i^g|Affree(ame)P 


167rm|mg 


the cross section simplifies to 


(Tg 


dfq 




q 


2m, 


-qvcosOqy^ X |Fdm(9)|^|/i^2(9)|^ • 


(A.7) 

(A.8) 

(A.9) 

(A. 10) 
(A. 11) 

(A. 12) 


A.2 Average rate in a dark matter halo 

The actual rate of excitation events, for a given transition and a given target electron, is found by 
multiplying Eq. (A. 12) by the DM density and averaging over the DM velocity distribution gx{v). 


Ri^ 2 = — d^v gx{v) avi^2 ■ 


m 


(A. 13) 
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In general, both the electron wavefunctions and the DM velocity distribution will not be spheri¬ 
cally symmetric. As noted in [9], the rate will then depend on the orientation of the target with respect 
to the galaxy. Here we ignore this interesting complication, and approximate the velocity distribu¬ 
tion as being spherically symmetric. We can then use the d^v integral to eliminate the (5-function in 
Eq. (A. 12), giving 

Ri^ 2 = Px/mx^ j ^ j ^ gx(^) Q(^ ~ ^min(g, AE;i^2)) X |Fdm(9)|Vi^2(9)P 

= Pxl^x^\- [ ■ (A.14) 

OTXPxe J q 

Here 7y(umin) has its usual definition, 

r]{v^\n) = j ~ dxi'^) - 

and Umin is a function of q and the energy transfer given by 

(XT? ^ ^^ 1^2 , 

q 

A.3 Ionizing an isolated atom 

For the purposes of connecting with previous work [9], in this subsection we consider ionization 
of electrons bound in isolated atomic potentials. We derive the ionization rate of such a system, 
assuming a spherical atomic potential and filled shells. This approximation was used in [9] to model 
a liquid xenon target material, and the results below reproduce Eqs. (5) and (6) of that paper. The 
full calculation of event rates in liquid xenon would require knowledge of electron wavefunctions in 
the dense, disordered xenon liquid. This is a more challenging calculation than for a semiconductor 
crystal, but can in principle be performed with similar methods - we leave this for future work. The 
corrections, however, can be argued to be small, lowering the ionization threshold and increasing the 
event rate. 

An electron ionized from an atom can be treated as being in one of a continuum of positive-energy 
bound states. These states are affected by the potential well of the atom, but can be approximated as 
free particle states at asymptotically large radii. We denote their wavefunctions as where 

I' and m! are angular quantum numbers, and k' is the momentum at asymptotically large radius. The 
energy of such a state is therefore Er = /2me- 

The ionization rate for such an atom is found by taking Eq. (A.14), summing over occupied 
electron shells, and integrating over the phase space of all possible ionized states. Since these are 
asymptotically free spherical-wave states, the phase space is 

■ ■ . . . V- fk'‘^dk' 1^ [ k'^dluEn 

ionized electron phase space = j = 2^ J ^ 

I'm' ' ' I'm' ' 

when the wavefunction normalization is, as in [9], taken to be 

~ ~ si/ 

{'4^k'l'm'\'4’klm} (2^) ^I'l^m'm ^2 k ) . (A.18) 
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Plugging this in, the ionization rate is given by 


R\, 


Px 

167r/r|g 


(A. 19) 


^ ^ J ^ Esi + k''^/2me))\FBM{q)\‘^\fi^k'i'm'{q)\^ ■ 

occupied I'm' ' ^ 

states 


where Esi is the binding energy of occupied state i. 

Since the potential is assumed to be spherically symmetric, and we are ionizing a full atomic 
shell, we can sum \fi^k'Vm'{^)\‘^ over initial and final angular momentum variables and the result 
cannot depend on the direction of q. This means we can define fhe dimensionless ionization form 
factor. 


\fion{k',q)f 


2k'^ 

(27r)^ 


E E 

occupied I'm' 
states 




(A.20) 


Affer applying fhis definition fo fhe previous equafion, we can replace fhe d?q infegral wifh Arrq^dq, 
giving 


dRion 
d In Efi 


Px/^x 


CTe 


qdq |-FDM(g)|^|/ion(fc',q)|^r/(nmin(Q',-Esi + jlrrie)) ■ 


(A.21) 


This reproduces fhe formulae given in [9]. 


A.4 Excitations in a semiconductor crystal 

In the periodic lattice of a semiconductor crystal, the electron energy levels form a complicated band 
structure, with an energy gap separating the filled valence bands and the unoccupied conduction bands 
(Fig. 4). Each possible electron level is labelled by a band index i and a wavevector k in the first 
Brillouin Zone (BZ). Due to the periodicity of the potential, the wavefunctions of these states are in 
Bloch form, 

f.^{x) = ^ ^ , (A.22) 

G 

where the G’s are the reciprocal lattice vectors. Here V is the volume of the crystal, and the wave- 
functions are taken to be unit-normalized, so that 

Y,\uiik + G)\^ = l (A.23) 

G 

(We use the relations f = (27r)^S^(k) and {27r)^6^{6) = V.) 
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With this form for the wavefunctions, the form factor Eq. (3.10) to excite from valence level {i k} 
to conduction level {i' k'} becomes 


= I E + G + G>,(t + G) ^ 


GG' 
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G' 


(27r)^d^(q — {k' + G' — k)) 
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u*,(k' + G + G')ui{k + G) 
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(A.24) 


(A.25) 


We define the term in the absolute square in Eq. (A.25) to be fr-r .,p 


[ik,i'k',G'y 


■f[ik,i'k',G'] 


'^uUk' + G + G')uiik + G). 


(A.26) 
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Inserting this into Eq. (A. 14), we can use the J-function to eliminate the d^q integral, giving 
_ Px 1 1 


R. 


Pie V ^ Q 


Yj ^ik)) I-^Dm(Q')|^ \f[ik,i'k',G']\ 


q=\k'+G'-k\ 

(A.27) 


The total excitation rate for an electron in level {i k} is found by summing Eq. (A.27) over all 
unfilled final energy levels i', 




^ Jbz (2vr)3 


(A.28) 


Note that we do not sum over final electron spins here as that sum has already been included in the 
definition of Ue- 

The total rate of excitation events in the crystal, i?crystab is given by summing Eq. (A.28) over 
all filled initial levels i, 

f Y(]3 

-^crystal = ‘^Yl (27r)3 ' (A.29) 

Here the extra factor of 2 is the sum over the two degenerate spin states of the filled valence bands. 
Putting this together gives the total excitation rate in a crystal. 
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crystal — 


^X Pie i^JBZ 


d^k d?k' 

(27r)® 


^^(^min(7, E.,p 
G' 




(A.30) 

where again q = \k' + G' — k\. Note that this is the total event rate for the whole crystal, and so it is 
appropriate that it is proportional to the volume V of the whole crystal. Since the dependence on the 
DM velocity distribution and interaction type are entirely encoded in rj and Fdm, which are functions 
only of the momentum transfer q and energy deposited E^, it is useful to insert delta-functions into 
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the above expression as follows: 

^xe 

f d^kd^k' 


crystal = / d In £"6 d In (vminCo', ^e)) i^DM (?)' 

^X ^xe j d. 




j^E,S{E, - E,J, + E.J) - |t' + G' - t|)|/[.j,„-, 3 /. (A.31) 

G' 


Using V = A'^ceiiUceii, where Uceii is the volume of the crystal’s unit cell and A^ceii is the number 
of cells, the differential rate can then be written in the form of Eq. (3.13), 

^ :^^celld^ea X ^ / dlnq (^r]{vrnm{q,Ee))^FBM{qf\fcryst!d{q,Ee)f , (A.32) 




n. 


xe 


where the crystal form-factor is defined as in Eq. (3.17), 

If ! p ^|2 _2vr2(«"leK;ell)"^ /■ Vce\\d?kVce\\d?k' 

|/crystal( 7 ,iiej| - ^27r)3 "" 


E-ME, - E.,j, + Ej) - It' - t + G'|)|/|,j,,s,_ 3 / . (A.33) 


G' 


B Derivation of inverse mean speed, r]{vmin) 


In this section, we will derive analytic expressions for r]{vmin)- For simplicity we assume that the DM 
velocity distribution, g-^{v^), takes the form of a Maxwell-Boltzmann distribution in the galactic rest 
frame, with a hard cutoff at the galactic escape velocity. In the Earth’s frame the velocity distribution 
then takes the form 

5x(F'x) = ^0 0(Uesc - \Vx + ve\) , (B.l) 

where is the DM velocity in the Earth frame, and ub is the Earth’s velocity in the galactic rest 
frame. We take vq = 230 km/s for the typical velocity, and Vesc = 600 km/s for the escape velocity. 
We take ve = 240 km/s for the average Earth velocity relative to the DM halo, adding (subtracting) 
15 km/s for the Earth velocity in June (December). 

The normalization factor K is determined by requiring f d^vg^iv) = 1, giving 


K = ngTr 





(B.2) 


Using these values, we obtain K = 6.75 x 10^^ [cm/s]^ or 2.50 x 10 ® in natural units. 
We then define the function ^^(umin). 


rjivmin) = / d^V^gxi^x) -©(^'x “ l^min) 

J Vx 

= ^ /'27rdcos6»c/uxUx ^ 


where cq = cos 6 is the angle between the velocity and the velocity of the Earth. We can explicitly 
solve Eq. (B.3), but need to consider two cases: 
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Figure 10. The minimum velocity, Vmin as a function of q for = 1 MeV, 10 MeV, and 1 GeV and = 1, 10, and 
40 eV. The grey shaded region indicates where Vmin > ve + Vesc and is thus prohibited. The vertical grey line indicates 
the maximum q value for a given E^ut = 70 Ry. 


1- Tmin < Vqsc '^E 

2- Tesc Ve < Tmin < ttesc “1“ Ve 


where Vesc-j Ve^ Vmin ^ 0* 

This gives us 
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2 / 2 

- Tmin + ve) + V^Vo 


Erf ( — I - Erf 

V ^0 


^^min — Ve 


Vo 


(B.5) 


where the subscript corresponds to the case number. Note that the two cases converge to the same 
value for Vmin = Vesc - VE- 


C Convergence of the Numerical Results 

In this section, we investigate the dependence of our calculation on the kinetic-energy cutoff, E'cut, 
(see Eq. (4.3)) and on the A:-point mesh. The choice of E^cut determines the maximum allowed three- 
momentum transfer q, which impacts the maximum Eg that we sample. Truncating the range of q 
can have more of an effect for high DM masses and high electron thresholds, since these two regimes 
depend on high values of q. On the other hand, the A:-points included in the mesh determine the 
computationally allowed values of q. A higher-resolution mesh is particularly important for low Ef. 
transitions to the bottom of the conduction band, and is therefore especially important for low DM 
masses and low electron thresholds. 


- 37 - 



















B 



[MeV] 


E, [eV] 



Q 


Figure 11. Left: Cross-section sensitivities for ionization thresholds of Qth = 5 and Qth — 11 electrons in silicon for 
Ecut = 30, 50, 70, 90, and 110 Ry (we take a mesh consisting of 27 fc-points). Note that most lines are on top of each 
other, demonstrating the weak dependence of (Te on Ecut- Right: Difference m the rate, 7?, for a given E^cut^ to that at 
Ecut = 70 Ry, Ro, in silicon for — 1 GeV. We see that the error in choosing Ecut = 70 Ry is < 0(1%) for the 
thresholds considered in this paper. 
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Figure 12. Left: Cross-section sensitivities for ionization thresholds of Qth = 5 and Qth — 11 electrons in germanium 
for Ecut = 30, 50, 70, 90, and 110 Ry (we take a mesh consisting of 27 fc-points). Note that most lines are on top of 
each other, demonstrating the weak dependence of CTe on Ecut- Right: Difference in the rate, R, for a given Ecut, to that 
at Ecut — 70 Ry, Ro, in germanium for = 1 GeV. The structure of the distributions arise from the effect of the 3d 
electrons. 


In Fig. 10, we show the dependence of Vmin in Eq- (3.12) on q for different and E^- The 
choice of Ecut (top axis) determines the range of q (bottom axis). In Fig. 11, we show the ratio of the 
rate for different values of Ecut to the rate at Ecut = 70 Ry for = 1 GeV for silicon targets. We 
see that the error in the rate with our choice of Ecut = 70 Ry is < 0(1%). The Ecut convergence 
in germanium is slightly worse due to the fact that we are solving for the 3d electrons instead of 
including them in the pseudopotential. This effect is greatest near the 3d shell energies (a few percent 
uncertainty), as seen in Fig. 12. In left plots of both Fig. 11 and 12, we see a step-like transition 
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Figure 13. Left: Cross-section sensitivity for Q = 1 and Q = 5 electrons in germanium for 27, 77, 137, and 243 fc-points. 
Right: The energy spectra using 27, 77, 137, and 243 fc-points for = 10 MeV and ae = 10“^^cm^. We see that the 
choice of the number of A:-points used in the mesh has an effect at low DM masses and low Q. 



Figure 14. The effects of our choice of A:-point mesh on our germanium results by perturbing the mesh with random 
shakes of amplitudes up to half the lattice-spacing. On the left, we look at the standard deviation of the shaken runs over 
the mean value as a function of DM mass. On the right, we look at the standard deviation of the shaken runs over the mean 
value as a function of for = 5 MeV. We see that our choice of A:-mesh spacing is accurate to a few tens of percent 
for masses above 1 MeV. 


between 30 Ry and the other curves for 1 le because 30 Ry is not a high enough energy cutoff to 
accurately calculate the rate for lie. The irregular behavior in the distributions on the right side of 
Fig. 12 are from the semicore electrons in Ge. We do not see the same behavior in Fig. 11, which 
considers Si and no semicore electrons. 

We investigate the effects of our choice of fc-point mesh on our results in two ways.^ First, we 
vary the number of fc-points in our mesh and find that there is sensitivity to our choice at low masses 
and thresholds, see Fig. 13. Second, we perturb each point on the mesh with a random shake of 
amplitude up to half the lattice-spacing so as to cover the entire fc-space. We use an energy cutoff 

*We do this only for the valence electrons without including the 3d-shell electrons, since energy of the latter are nearly 
constant as a function of k. In any case, the 3d-shell electrons are important at large E^, while the choice of A:-point mesh 
is important only at low E^. 
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of Ecut = 70 Ry and 243 A:-points. The amplitude of our perturbations is AA; = 0.08 a.u. as the 
lattice spacing for 243 fc-points is 0.17 a.u. We run 10 independent simulations and plot the results 
in Fig. 14. We find that our choice of A:-point mesh does not appreciably affect our results for masses 
above IMeV. 


D The importance of the 3d-shell in germanium 

The importance of the 3d-shell electrons in germanium are illustrated in Fig. 15. We see that they 
dominate the rate at high recoil energies and thus for high thresholds. We discuss this in more detail 
in S6.2. 


E, [eV] 

0.67 3.6 6.5 9.4 12.3 15.2 18.1 21. 23.9 26.8 29.7 32.6 35.5 38.4 41.3 




[MeV] 

Figure 15. The importance of the Sd-shell electrons in germanium. Left: Spectrum of events as function of the ionization 
signal Q, for Tdm = 1 and CTj. = 10~^^cm^. The thick, upper lines show the rates including the 3d-shell electrons, 
while the thinner, lower lines include only the valence electrons (the thick and thin lines overlap for Q < 8). The shading 
highlights the difference between the two. Right: The cross-section reach in germanium with a 1 kg-year background-free 
exposure, with an ionization threshold of Qth = 5 for the dashed, lower lines and Qth = 10 for the solid, upper lines. 
The lower (upper) line of the shaded region show results with (without) the 3d-shell electrons (these overlap for the 5e 
threshold). 


E A Monte Carlo model of secondary scattering 

In the main results of this paper, we modeled the ionization response of a target crystal with the linear 
treatment described in §5. For comparison, here we attempt to mock-up the secondary scattering with 
a Monte Carlo model, following [101]. The deposited energy Eg, is randomly split between an initial 
electron and hole. In each following step, each electron or hole with energy above a threshold E'ion 
then generates an extra electron-hole pair, with the energy being randomly split between the three 
particles. This is iterated until all particles have energy less than Sion- 

The random energy splittings follows a distribution that weights all phase space volume equally, 
with the density of states assumed to grow as y/S above and below the bandgap, as in a simple 2- 
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Figure 16. Probability distribution of the ionization signal Q for a given energy deposition Es in germanium. The left 
plot shows the distrihution of Q for the indicated fixed Ee, while the right plot shows the probability to get a given Q with 
varying E^. Solid, filled lines show the cascade model discussed in §E, while dashed lines show the linear model described 
in §5 and used for our main results. 
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Figure 17. Fractional increase in the rate when modeling the secondary scattering with the cascade Monte Carlo instead of 
the linear model. (Explicitly, the y-axis is (i?MC — Knaive) / Rna.ive-) The rate here is the rate of event passing the ionization 
threshold Qthr, for a DM mass of 1 GeV. See §E for more details. 


band free electron/hole system. Explicitly, for the initial 1^2 splitting the probability distribution 
for energy Eq to split into energies Ei and E2 has the form 

dP (X El y/ E2 S{Eq — El — E2 — -Egap) dEi dE 2 , 

while for the subsequent 1 —)• 3 splittings it has the form 

dP oc ^/Ei^/E 2 -\/e^^{Eq — El — E2 — E^ — E'gap) dEi dE2 dE^ , 

where electon/hole energies are measured above/below the upper/lower edge of the band gap. We 
ignore phonon losses during the cascade - these are understood to be quantitatively fairly small, and 
should not affect the qualitative conclusions. 

The output of the Monte Carlo model is a probability distribution P{Q\Ee) to get ionization Q 
given a deposited energy Eg. Given the band-gap energy of E^gap = 0.67 eV (1.11 eV) in germanium 
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(silicon), we find that E^on = 2.67 eV (3.1 eV) reproduces the measured values of e for high energy 
recoils (see Eq. (5.2)). The distributions for both elements have Fano factors of F « 0.1 for all 
energies above ~ 10 eV, consistent with measurements. We illustrate the probability distributions in 
Fig. 16. In Fig. 17, we show the effect on the event rate of using this model rather than the naive 
linear model of §5. For thresholds of 2 to 4 electron-hole pairs, downward fluctuations reduce the 
rate compared to the naive estimate. For higher thresholds, occasional upward fluctuations combined 
with the steeply falling recoil spectrum lead to an increase in the rate. However, the two models are 
consistent within a few tens of percent. 


F Review of Density Functional Theory and Pseudopotentials 

In this appendix, we review the formalism of density functional theory (DFT), explain in more detail 
the approximations used in the computation of the wavefunctions, and further explain the numerical 
methods. 


F.l Electronic structure within DFT 

Non-relativistic electrons interacting electrostatically with fixed nuclei are described by the electronic 
structure Schrddinger equation 


1 

2mp 


EvJ-E 


Zje^ 


aJ 


-R, 


+ E 

a<p 


-F3| 


T'i (fi,... ,r)v) = ejfl'i (fi,... ,fAr), (F.l) 


where a,j3 = 1, 2, • • • , label electrons, I = 1, 2, • • • , M labels nuclei, and Z/ is the atomic 
number of nuclei I. The first term in the Hamiltonian is the electron kinetic energy T, the second term 
is the Coulomb electron-nucleus attraction 14xt and the third term is the electron-electron Coulomb 
repulsion Vee- The constant nuclei-nuclei term has been omitted. Even though the question is well- 
posed, obtaining the many-electron wavefunctions (fq,..., Ttv) computationally is an extremely 
difficult task because of the exponential scaling of the problem with the number of electrons N. 
This method becomes then helpless for applications of interest, so in practice one needs to resort to 
approximate methods. 

DFT is a reformulation of the interacting quantum many-body problem in terms of functionals 
of the particle density n(r). For the case of electrons, the Hohenberg-Kohn theorems [91] imply 
that all properties of the interacting system are determined once the ground state electron density 
is known. Minimizing an energy functional E [n] will provide the ground state density no(r) and 
the ground state energy Eq. Unfortunately this energy functional is not known in general. The 
Kohn-Sham method [92] overcomes this obstacle by replacing the description strictly in terms of 
functionals for a wavefunction formulation: the system of interacting electrons with Hamiltonian 
H = T + Vee + Kxt is mapped into a system of independent electrons under the presence of an aux¬ 
iliary potential H = T + Vaux + Kxt which produces the same ground state density as H. This is of 
great advantage because, once this mapping is built, one has to solve the much simpler independent- 
particle system in order to obtain Eq and no(r). However, this comes at the expense of having to use 
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an approximate auxiliary potential 14ux- Typically, 14ux is split into the mean-field Hartree potential 
i4tartree(i^ = ^ f d^'f^n{f*)/ |f — f^| and the so-called exchange-correlation potential 14 c, where the 
approximations are imposed. Once an approximation for 14c has been chosen, the non-interacting 
electron Schrodinger equation 


2m, 


+ 14xt(T) -|- 14[artree(?*) + 14c(f0 




(F.2) 


which are known as the Kohn-Sham equations, are solved to get the auxiliary Kohn-Sham wavefunc- 
tions 'ipi{r). From these, the density can be obtained as n{f) = Ylifi where fi are the 

occupation numbers {fi = 2 for spin-unpolarized systems) as well as the ground state energy by eval¬ 
uating the energy density functional^ E[n] = T [n] -|-F^ext [^] + -^Hartree [n] + Exc as well as a set 

of wavefunctions ifi and eigenenergies Si. This problem is solved self-consistently until convergence 
is reached. 

Expanding the wavefunctions in a finite plane-wave basis with elements labeled by the vectors G 
and G', the Kohn-Sham equations become the matrix equations 


where the Hamiltonian is 


''^HgQ,{k)ui{k + G) = Ei{k)ui{k + G). 

G' 


k + G 


H 


k + G') = 


2m e 


k + G 


6.., + ViG-G') 


(F.3) 


(F.4) 


It should be noted that, since the potential is local, its reciprocal space form does not depend on k. 
Furthermore, the Kohn-Sham equations in reciprocal space Eq. (F.3) decouple different fc’s, so the 
eigenvalue problem can be carried out independently at each k. 

Despite all the successes of DFT, several notable shortcomings are known today. The most rele¬ 
vant one for us is that DFT is known to give an incorrect band gap. This is due to a discontinuity in the 
DFT exchange-correlation potential dExc/dn{f) when electrons are added above the gap [115, 116]. 
There are methods based on many-body perturbation theory to improve the DFT band gap and band 
shapes, such as the GW method [117]. However, since the largest contribution to the scattering rate 
comes from the low energy excitations, we choose to follow an empirical “scissor correction” ap¬ 
proach [95, 96]. In this approach a rigid shift is imposed on the conduction bands with respect to 
the valence bands in order to set the band gap to the experimental values of 1.11 eV for silicon and 
0.67 eV for germanium [97]. It is worth noting that the semiconductor band gap features a tempera¬ 
ture variation of around 10 meV [94], but we are choosing the room temperature band gap values for 
our calculation. 


F.2 Energy Density Functionals 

In order to be able to use DFT, a choice for the exchange-correlation functional Exc [n] is required. 
The Local Density Approximation (LDA) [92] has been remarkably successful because of its sim¬ 
plicity and transferability. In this method the exchange-correlation energy functional is based only 

®The connection between an energy functional and its corresponding local potential is E\n\ = J (Efn{r) V(r). 

'®The kinetic energy is calculated from the Kohn-Sham wavefunctions as T = ll2ms \ W%j}i\^. 
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on physical considerations and is approximated locally by the energy of a homogeneous electron gas 
(HEG) with the following density: 

N = J + 6™ ^(n(r-l)] . (F.5) 

The HEG exchange [118] and correlation [119] energy functionals are well established. There are 
some faults in the EDA which are known to be most dramatic where the electrons are highly local¬ 
ized and exchange repulsions are significant. In order to correct for that, the Generalized Gradient 
Approximations (GGA) introduce a dependence on the density gradient in the exchange-correlation 
energy density 

[n] = J d^rn{^ e°f'^(n(r), |Vn(r)|). (F.6) 

In this work we choose the well-established PBE functional [93] which is known to produce a broad 
set of properties of materials to accuracies of order a few percent [120]. Since EDA functionals tend 
to underestimate the energies of excited states compared to GGA functionals, we find a difference in 
cross-secfion sensifivify of around 10-20%, wifh a larger difference af higher fhresholds. 

F.3 Pseudopotentials 

The valence electrons are responsible for the formation of interatomic bonds and their wavefunctions 
are in general delocalized, spanning over interatomic distances. The core electron wavefunctions, 
however, are very localized around the nucleus and they barely change from the isolated atom to 
the condensed matter phase. This fact allows to use the atomic core electron wavefunctions in the 
condensed matter phase by replacing the bare positive nuclear Coulomb potential and the negative 
Coulomb potential generated by the core wavefunctions with a pseudopotential in the Kohn-Sham 
problem Eq. (E.2). The advantage is two-fold: first, the number of electrons in the problem is reduced 
to the number of valence electrons and second, the only wavefunctions to be calculated are valence 
wavefunctions which, since they are rather smooth, do not require as fine a grid to represent them as 
a core electron wavefunction would, thus improving the computational efficiency. In this work we 
use Vanderbilt-type ultrasoft pseudopotentials [94]. The pseudopotential for Si includes the 3s and 3p 
electrons in the valence, while in the case of germanium, we use a pseudopotential which includes the 
3d, 4s and 4p electrons in the valence. 
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